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Background 

Ordnance  detection  and  discrimination  techniques  are  typically  based  on  two- 
dimensional  representations  of  geophysical  data  such  as  contoured  grids.  These  grids  are 
created  by  interpolating  data  collected  along  transect  lines.  Without  the  use  of  wheeled 
carts  (which  pose  their  own  set  of  unique  restrictions)  the  assumption  that  all  data  are 
collected  at  a  common  height  above  the  ground  is  routinely  violated.  Differences  in 
sensor  height  between  adjacent  lines  cause  systematic  and  pseudo-random  distortions  in 
the  gridded  data  which  are  not  based  on  geophysical  changes.  The  resulting  anomalies 
are  easily  misinterpreted,  or  require  filters  which  reduce  the  reliability  of  inversion 
routines.  Field  efforts  to  collect  data  of  the  quality  required  for  adequate  inversion  are 
painstakingly  conducted  to  achieve  a  measure  of  unifonnity.  The  time,  effort,  and 
training  required  by  these  efforts  have  seldom  been  translated  into  cost-effective  project 
implementation. 


Objective 

The  objective  of  this  research  is  to  determine  some  of  the  benefits  accrued  by  collecting 
and  processing  three-dimensionally  positioned  data  (3D)  for  unexploded  ordnance 
(UXO)  detection  and  discrimination  versus  traditional  two-dimensional  (2D)  approaches. 


Technical  approach 

The  technical  approach  taken  here  was  strictly  model-based  in  that  no  geophysical  data 
were  collected  at  this  stage  of  the  project.  Although  the  techniques  are  equally  applicable 
to  electromagnetic  as  well  as  magnetic  data  sets,  only  magnetic  data  were  modeled  at  this 
stage.  2D  and  3D  collection  scenarios  were  devised  and  tested  over  a  variety  of  single 
and  double  target  models.  Simulated  targets  included  60mm  and  20mm  shells  at 
relatively  deep  burial  depths  and  close  proximity.  Moderate  levels  of  noise  were  added 
to  the  positions  and  to  the  synthetic  magnetic  readings. 

In  this  project,  3D  processing,  visualization,  detection  and  discrimination  techniques 
were  developed  as  analogs  to  the  2D  techniques  most  commonly  in  current  use. 
Comparisons  were  made  between  the  2D  and  3D  techniques  to  evaluate  the  relative 
merits  of  the  3D  approach. 

A  new  collection  procedure  was  devised  to  maximize  the  utility  of  the  3D  approach.  This 
is  principally  a  “reacquisition”  or  “cued  investigation”  approach  in  which  an  operator  is 
directed  to  an  anomalous  location  on  a  dig  list  and  proceeds  to  collect  data  over  the  target 
area  by  waving  the  sensor  randomly  over  the  center  point.  This  is  similar  to  the  approach 
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currently  used  by  OE  Technicians  with  analog  magnetometers,  only  in  this  case  the 
instruments  would  be  digital  and  capable  of  recording  3D  positioned  data.  Eventually, 
the  processing  and  inversion  techniques  developed  here  can  be  implemented  on  a  real¬ 
time  basis  in  order  to  minimize  field  time. 

Processing  techniques  were  extended  to  3D  using  commercial  software.  A  variety  of 
“gridding”  routines  were  tested  for  their  ability  to  reliably  interpolate  potential  field  data 
and  an  optimal  technique  was  chosen.  Hanning  filter  coefficients  were  also  calculated 
and  programmed  to  smooth  the  data  in  three-dimensions. 

Visualization  techniques  are  irrelevant  to  numerical  inversion  routines,  but  are  required 
for  operator  input.  For  example,  without  a  simple  visualization  technique  the  operator 
cannot  tell  whether  to  invert  for  one  or  more  objects.  A  variety  of  standard  visualization 
techniques  were  examined,  but  the  most  effective  was  the  Gradient  String.  This  is  a 
unique  concept  devised  specifically  for  this  project.  A  Gradient  Strings  traces  a  line  of 
maximum  magnetic  gradient  through  the  three-dimensional  data  space  above  a  target. 

This  reduces  the  solid  cube  of  data  into  a  single  strand  which  points  directly  to  the  target 
or  targets. 

Discrimination  techniques  included  a  variety  of  3D  analogs  to  2D  techniques.  For  the 
purposes  of  this  report,  discrimination  is  differentiated  from  detection  in  that  detection  is 
limited  to  identifying  the  number  of  targets,  whereas  discrimination  attempts  to  determine 
specific  parameters  such  as  depth,  size  and  orientation.  In  particular,  higher  order 
gradients  were  used  to  estimate  depths  (Salem  et  al,  2002),  as  was  a  simple  Euler 
deconvolution.  The  MTADS-DAS  dipole  model  code  was  used  to  estimate  depth  and 
size.  The  commercial  ModelVisionPro  software  was  used  to  estimate  depth,  size  and 
orientation.  Discrimination  in  this  study  does  not  include  differentiating  UXO  from  non- 
UXO. 


Summary 

The  3D  acquisition  and  processing  methodology  was  shown  to  be  better  than  the  standard 
2D  method  for  detection  and  location  of  closely  spaced  targets  in  all  cases.  The  Gradient 
String  approach  detected  all  of  the  targets  and  located  them  in  (x,  y)  more  accurately  than 
the  comparable  technique  of  2D  analytic  signal  peaks.  Depth  estimates  using  the  analytic 
signal  ratio  and  Euler  deconvolution  were  five  times  more  accurate  using  the  3D  data. 

The  most  accurate  technique  was  the  3D  ModelVision  spherical  inversion,  which 
averaged  location  errors  of  9cm  and  depth  errors  of  2cm. 

Tests  showed  that  detection  and  inversion  of  2D  noise-free  data  was  comparable  to  the 
3D  noise-free  data  in  most  cases.  This  is  not  unexpected  for  potential  field  data.  The 
introduction  of  noise,  however  (1cm  of  positional  noise  and  O.lnT  of  magnetometer 
noise),  created  significant  problems  for  the  2D  data.  In  general,  the  data  were  insufficient 
to  invert  or  (in  some  cases)  detect  multiple  objects  even  with  tight  line  spacing.  The  3D 
data  also  suffered  from  the  noise,  particularly  at  low  heights  where  gradients  were 
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strongest,  but  detection  and  location  capabilities  were  almost  always  2-5  times  better  than 
the  comparable  2D  data. 

The  Gradient  String  technique  proved  to  be  an  excellent  method  of  detecting  and 
visualizing  multiple  targets  within  the  data  space.  The  various  inversion  techniques 
generally  improved  the  2D  results,  but  in  almost  all  cases  the  3D  processing 
outperformed  the  comparable  2D  processing.  The  3D  data  provided  enough  additional 
information  to  improve  detection  and  depth  estimates  by  up  to  5  times,  but  not  enough  to 
determine  more  complex  parameters  such  as  dimensions  and  orientation. 


Table  1:  Average  location  errors  by  method.  Failures  in  the  dipole  inversion  are  generally  the  result 
of  insufficient  signal  strength. 


2D  data  set 

3D  data  set 

Peak  detection 

horiz-error 

0.18  m 

0.09  m 

AS  ratio  depth 

vert-error* 

0.21  m 

0.10  m 

Euler  depth 

vert-error* 

0.24  m 

0.06  m 

Dipole  inv 

radial-error 

0.12  m  (9  failures) 

0.28  m  (1  failure) 

Sphere  inv 

radial-error 

0.14  m 

0.09  m 

Ellipse  inv 

radial-error 

0.17  m 

0.18  m 

*  vertical  error  listed  1 

lere  is  at  location  with  ’ 

horizontal  error  given  for  peak  detection 

(i.e.  Euler  2D  radial  error  =  V(  0.182  +  0.247)  =  0.30m) 

Several  lines  of  development  are  suggested  by  the  results  of  this  project.  Implementation 
of  feedback  mechanisms  and  an  extension  of  the  3D  techniques  to  electromagnetic 
methods  are  recommended  to  extract  the  additional  features  such  as  orientation  from  the 
data.  More  advanced  data  segmentation  techniques  are  also  recommended  in  order  to 
more  effectively  resolve  model  parameters.  Discrepancies  between  the  various  software 
packages  for  ellipsoid  model  signatures  were  found  and  need  to  be  resolved.  Test-stand 
measurements  to  verify  signature  shapes  and  field  trials  to  verify  performance  standards 
are  required.  Similarly,  tests  to  improve  inversion  results  for  target  shape  and  orientation 
and  for  discrimination  between  UXO  and  non-UXO  are  still  required  to  fully  implement 
this  3D  technique. 
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Accomplishments 

Introduction 

This  project  is  broken  down  into  several  parts,  each  with  its  own  method  and  conclusions, 
and  most  flowing  sequentially  from  one  to  the  other.  For  reporting  continuity,  these 
topics  are  covered  individually  rather  than  providing  method  for  all  parts  and  then  results 
for  all  parts.  These  topics  include: 

•  Project  Conception 

•  Data  Collection  and  Modeling 

•  Gridding,  Filtering  and  Processing 

•  Visualization  and  Gradient  Strings 

•  Detection 

•  Depth  Estimation 

•  Euler  Deconvolution 

•  Dipole  Inversion 

•  Solid-body  Inversion  (Sphere) 

•  Solid-body  Inversion  (Ellipsoid) 
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Project  Conception 

Many  geophysical  data  sets,  collection  procedures  and  presentation  techniques  are 
described  as  3-dimensional.  For  the  purposes  of  clarity,  a  distinction  must  be  made 
between  the  number  of  dimensions  of  the  positioning  data  recorded  during  the  survey  and 
the  number  of  dimensions  of  the  presentation  fonnat.  The  traditional  survey  records  two 
dimensions  of  position  data  (x,  y)  and  presents  the  geophysical  results  as: 

•  ID  profiles  (a  projection  against  a  single  axis  of  time  or  horizontal  position,  with 
the  data  as  the  second  axis), 

•  2D  grids/contours  (a  projection  against  two  horizontal  position  axes  using 
contours  or  colors  to  represent  the  data  axis), 

•  or  3D  surfaces/curves/nets  (a  projection  against  two  horizontal  position  axes, 
using  the  illusion  of  a  third  vertical  axis  in  data  units  to  represent  the  geophysics, 
often  including  color  and/or  contours  for  emphasis  of  the  data). 

Even  with  the  recording  of  3D  positional  data  (x,  y,  height),  none  of  the  aforementioned 
presentation  techniques  display  or  accommodate  the  addition  of  height  infonnation. 
Generally,  the  height  is  assumed  to  be  constant  and  the  data  are  forced  onto  a  2D  grid. 

We  must  also  recognize  the  limitations  of  all  two-dimensional  media  (paper  and  screens). 
Even  3D  surfaces  are  technically  2D  presentations  using  illusions  of  perspective  and 
shading  to  appear  3D.  It  is  extremely  difficult  to  represent  3D  points  on  a  2D  screen  or 
map,  and  this  does  not  yet  include  a  dimension  for  the  geophysical  data.  Occasionally, 
the  perspective  of  a  3D  presentation  is  rotated,  thus  adding  time  to  enhance  the  illusion  of 
dimension.  Exchanging  our  2D  media  for  3D  holographic  presentations  would  still 
require  some  form  of  transparency  plus  rotation  of  the  image  or  movement  of  the 
observer  in  order  to  fully  interpret  it,  just  as  any  solid  object  would.  A  3-dimensional 
representation  therefore  presents  several  fundamental  obstacles  to  the  observer. 

Multi-dimensional  data  sets  require  that  we  either  reduce  the  information  to  its  most 
relevant  fonn  for  presentation,  or  think  and  present  data  in  terms  of  matrices  and 
functions  with  n-dependent  variables.  This  applies  equally  to  3D  positioned  data  from  a 
single  sensor,  multiple  responses  from  a  single  sensor  (ie  EM  time  gates),  or  integrated 
multi-sensor  data  sets.  Although  this  makes  them  ideal  candidates  for  computer  analysis, 
which  is  based  almost  entirely  on  linear  matrices,  the  human  factor  is  still  required  for 
control  and  judgment,  hence  the  requirement  for  visualization. 

A  single  geophysical  measurement  at  one  point  in  space  rarely  has  any  significance.  A 
measurement  requires  spatial  correlation  with  neighboring  points  (ID  profile)  to  give 
data  meaning.  Similarly,  consecutive  profiles  allow  the  data  to  be  spatially  correlated  in 
two  dimensions  -  down  line  and  across  line.  The  traditional  approach  to  geophysical 
mapping  therefore  is  to  collect  data  along  parallel  transects  (ID)  at  a  nominal  height  and 
present  the  results  as  a  grid  (2D).  Each  added  dimension  brings  a  new  level  of  relevance 
to  the  interpretation,  but  there  are  drawbacks. 
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Single  profiles  present  the  data  with  great  resolution  and  fidelity  since  they  are  not 
subject  to  gridding  interpolation,  but  can  be  difficult  to  interpret.  The  most  obvious 
shortfall  is  the  lack  of  information  to  the  sides  of  a  profile  anomaly.  By  presenting  the 
data  with  reference  to  neighboring  points  on  both  sides,  the  2D  perspective  makes  it 
much  simpler  to  discriminate  discrete  targets  from  broader  anomalies.  This 
simplification  is  an  important  part  of  the  interpretation  process. 

Data  processing  techniques  based  on  2D  data  sets  also  have  significant  advantages  over 
their  ID  counterparts.  There  are  many  advantages  to  analyzing  data  that  have  been 
spatially  correlated  in  two  dimensions.  Spatial  filters,  for  example,  can  remove 
uncorrelated  noise  thereby  enhancing  the  signal-noise  ratio.  FFT  processes,  such  as 
analytic  signal,  can  also  enhance  features  more  reliably  and  with  fewer  assumptions  when 
applied  to  grids  than  to  profiles. 

The  disadvantage  of  2D  presentations  is  that  they  assume  that  all  of  the  data  were 
collected  on  a  common  plane,  although  this  is  rarely  true.  Whenever  a  data  point  is 
collected  at  different  height  than  the  adjacent  data  points,  the  gridding  process  creates  an 
artifact  in  the  grid.  Depending  on  the  severity  and  distribution  of  these  artifacts,  any 
interpretation  based  on  the  gridded  data  will  be  more  or  less  flawed.  As  a  result,  many 
additional  2D  data  processing  steps  must  be  incorporated  into  the  work  flow  in  order  to 
correct  for  the  “common  plane”  assumption  and  produce  a  viable  interpretation.  In 
addition  to  the  gridding  process,  which  generally  desamples  data  down-line  and 
interpolates  data  between  lines,  filtering  or  deletion  of  data  points  must  be  applied  to 
produce  a  satisfactory  grid.  In  some  cases,  these  additional  steps  may  reduce  data  quality 
more  than  they  enhance  it. 

The  impact  of  variable  sensor  height  on  gridded  data  has  been  addressed  by  many  others 
and  the  effects  on  both  magnetic  and  electromagnetic  data  are  well  documented. 
Numerous  studies  have  been  done  to  limit,  quantify  or  correct  these  effects.  Techniques 
include  wheeled  carts  to  control  height,  altimeters  and  other  positioning  systems  to 
measure  height,  and  various  upward  or  downward  continuation  algorithms  to  correct 
readings  to  a  nominal  height.  Advances  in  positioning  systems  and  computing 
capabilities  now  allow  collection  and  storage  of  geophysical,  position  and  height  data  in  a 
single  database.  Software  packages  for  data  inversion  are  taking  advantage  of  increased 
power  and  are  inverting  data  with  height  information.  Although  the  technology  has 
incorporated  several  aspects  of  3D  data,  all  of  these  are  attempts  to  solve  the  problem  of 
continually  forcing  3D  data  into  a  2D  format. 

To  date,  the  best  inversion  solutions  ignore  the  2D  format  altogether  and  work  entirely 
with  profile  data.  Although  this  work  flow  avoids  many  of  the  problems,  it  does  not  take 
advantage  of  any  of  the  positive  elements  of  2D  processing  such  as  noise  reduction  or 
FFT  algorithms.  Processing  benefits  are  limited  to  one  dimension  of  spatial  correlation 
(down-line).  Only  filters  and  processes  that  can  be  applied  to  profiles  are  used. 
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The  only  way  to  take  advantage  of  spatial  correlation  of  data  collected  with  3D 
positioning  is  to  process  it  in  three  dimensions.  Not  only  does  this  effectively  deal  with 
the  problems  associated  with  changes  in  height,  but  it  offers  the  potential  of  spatially 
correlating  the  data  vertically  as  well  as  horizontally.  This  also  provides  an  opportunity 
to  enhance  the  state-of-the-art  inversion  techniques  by  incorporating  the  benefits  of 
spatial  correlation  prior  to  inversion. 

This  project  was  therefore  designed  to  examine  the  steps  necessary  to  maximize  the 
utility  of  3D  positioned  data,  to  develop  the  groundwork  for  some  of  the  tools  that  may 
be  required,  and  to  quantify  the  benefits  accrued  through  such  a  process. 
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Data  Collection  and  Modeling 

For  logistical  and  testing  purposes,  this  project  assumed  that  the  data  were  collected  as 
part  of  a  cued  investigation  (to  follow-up  a  previously  detected  anomaly)  using 
magnetometer  sensors.  The  techniques  are,  however,  equally  applicable  to  traditional 
survey  collection  and  to  electromagnetic  data  after  additional  experimentation  with  3D 
gridding  and  inversion  techniques.  After  an  initial  geophysical  survey,  specific 
anomalies  are  identified  for  additional  follow-up.  Prior  to  excavation,  the  operator 
collects  a  “cloud”  of  data  points  at  irregularly  spaced  intervals  around  the  pre-determined 
location.  The  “cloud”  of  data  is  then  interpolated  into  a  3D  lattice  for  processing  and 
analysis  (Figure  1). 

Forward  models  were  created  to  develop  and  test  new  3D  algorithms  and  to  test  the 
comparison  between  2D  and  3D  capabilities.  Initial  testing  of  new  routines  and  software 
was  done  with  simple  dipole  models.  For  the  final  evaluation  of  the  3D  techniques,  new 
models  were  created  using  the  MAGMOD  code  as  the  basic  kernel.  It  had  the  advantage 
of  modeling  ordnance-like  sources,  and  was  independent  of  the  various  inversion  codes 
being  tested.  The  user  interface  and  output  fonnats  of  the  original  Matlab  code  were 
revised  to  accommodate  the  specific  project  requirements.  This  included  acceptance  of  a 
list  of  3D  measurement  points,  output  of  results  in  3D,  grid  orientations  other  than 
magnetic  North,  and  modeling  multiple  targets. 

The  final  list  of  models  included  single  and  closely  spaced  ordnance  at  various  depths 
and  orientations.  The  various  target  parameters  are  listed  below.  Non-varying 
parameters  are  listed  in  Table  2,  variable  parameters  in  Table  3.  Scale  representations  of 
the  models  against  the  2D  color  contours  are  shown  in  Appendix  B-l. 

Target  locations  were  chosen  to  represent  a  scenario  in  which  a  single  anomaly  may  be 
reasonably  detected  on  a  standard  survey.  Line  spacing  and  noise  levels  were  set  to 
simulate  a  search  for  small  UXO  (20mm)  on  a  mixed-use  range.  The  results  of  the  2D 
simulation  show  an  anomaly  that  would  be  selected  for  further  follow-up.  The 
reacquisition  phase  should  be  sufficient  to  determine  if  there  is  more  than  one  target 
contributing  to  the  original  anomaly  and  its/their  basic  parameters  (position,  depth,  size). 
Targets  were  placed  close  enough  together  to  overlap  signatures  at  typical  sensor  heights, 
but  still  have  distinct  positive  peaks  that  would  be  considered  separate  hits  with  a  0.5m 
search  radius.  Burial  depths  were  guided  by  maximum  penetration  limits  estimated  by 
USAESCH.  The  deepest  60mm  targets  were  placed  at  90%  of  the  maximum  penetration 
depth,  while  the  20mm  targets  were  placed  at  approximately  25%  of  maximum  depth. 
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Table  2:  Non-varying  target  parameters  for  final  set  of  models.  Target  dimensions  were  obtained 
from  ordnance  tables.  Relative  permeability  is  the  default  used  by  the  Magmod  code.  Remanence 
has  been  set  to  zero,  and  the  effects  of  demagnetization  and  octupole  components  have  been  included. 


60mm  equivalent  sphere 

Diameter 

0.102m 

Relative  permeability 

500 

20mm  equivalent  sphere 

Diameter 

0.034m 

Relative  permeability 

500 

60mm  ellipsoid 

Major  axis  length 

0.200m 

Semi-major  axis  length 

0.060m 

Minor  axis  length 

0.060m 

Relative  permeability 

500 

20mm  ellipsoid 

Major  axis  length 

0.069m 

Semi-major  axis  length 

0.020m 

Minor  axis  length 

0.020m 

Relative  permeability 

500 

Earth’s  field 

Strength 

50,000nT 

Inclination 

70°N 

Declination 

0°E 

Table  3:  Variable  target  parameters  for  final  set  of  models.  Negative  depth  is  below  ground. 
Azimuth  is  the  horizontal  angle  between  the  major  axis  and  grid  North  (positive  CW).  Dip  is  the 
vertical  angle  between  the  major  axis  and  horizontal  (positive  down). 


Model 

Type 

X 

Y 

Depth 

Azirn 

Dip 

1 

60mm  sph 

-0.1 

-0.2 

-0.2 

0 

0 

2 

60mm  sph 

-0.1 

-0.2 

-0.6 

0 

0 

3 

60mm  sph 

-0.1 

-0.2 

-0.6 

0 

0 

20mm  sph 

+0.1 

+0.3 

-0.05 

0 

0 

4 

60mm  ellip 

-0.1 

-0.2 

-0.2 

0 

0 

5 

60mm  ellip 

-0.1 

-0.2 

-0.6 

0 

0 

6 

60mm  ellip 

-0.1 

-0.2 

-0.6 

0 

90 

7 

60mm  ellip 

-0.1 

-0.2 

-0.6 

0 

0 

20mm  ellip 

+0.1 

+0.3 

-0.05 

0 

0 

8 

60mm  ellip 

-0.1 

-0.2 

-0.6 

0 

90 

20mm  ellip 

+0.1 

+0.3 

-0.05 

0 

90 

9 

60mm  ellip 

-0.1 

-0.2 

-0.6 

110 

65 

20mm  ellip 

+0.1 

+0.3 

-0.05 

45 

40 
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Data  sampling  for  the  2D  and  3D  data  sets  were  based  on  traverse  lines  and  cued 
investigations  respectively.  For  the  2D  set,  lines  were  spaced  0.5m  apart  at  a  nominal 
0.35m  height.  Five  lines  were  sampled,  covering  +/-  lm  from  the  target.  Lines  1,3,5 
were  set  exactly  at  0.35m  (including  the  line  closest  to  the  targets),  while  Line  2  was  set 
at  0.40m  and  Line  4  was  set  at  0.30m.  This  reflected  the  inherent  irregularity  of  adjacent 
line  heights  without  unduly  biasing  the  overall  average  or  the  line  closest  to  the  targets. 
Down  line  data  spacing  was  set  at  0.05m,  which  corresponded  to  a  slow  walking  speed  of 
0.5m/s  at  a  10Hz  sample  rate.  Although  roughly  half  nonnal  walking  speed  (and  twice 
the  data  density),  this  was  consistent  with  the  level  of  care  required  to  maintain  a  0.5m 
line  spacing. 

The  3D  set  was  based  on  a  0.05m  sample  increment  over  a  range  of  +/-  lm  horizontally 
from  the  center  point  (0,  0)  and  from  0. lm-0.6m  height.  A  random  selection  of  300  data 
points  were  taken  from  this  set  with  a  bias  toward  those  closest  to  (0,  0,  0).  This 
simulated  a  30  second  collection  period  at  10Hz  sample  rate  over  the  target  of  interest. 

Model  results  for  the  2D  and  3D  point  sets  were  generated  using  the  modified 
MAGMOD  code.  Normally  distributed  noise  of  0.01m  was  added  to  the  positions  in  (x, 
y,  z)  and  0.1  nT  noise  was  added  to  the  magnetic  field.  No  additional  noise  was  added  to 
simulate  heading,  diurnal  or  operator  error.  These  effects  would  only  impact  the  2D  data, 
since  the  3D  data  would  have  a  constant  heading,  the  diurnal  would  vary  only  slightly 
over  the  30sec  collection  period,  and  the  operator  would  remain  stationary. 

In  total,  there  were  205  data  points  in  the  2D  data  sets  and  300  points  in  the  3D  data  set. 
Although  this  gave  a  slight  advantage  to  the  3D  data  set  in  terms  of  the  number  of  raw 
data  points,  it  was  the  premise  of  this  research  that  the  collection  methodology  provides 
the  greatest  advantage.  These  are  the  data  sets  used  for  the  inversion  comparison.  Figure 
2  illustrates  the  basic  data  distribution  for  the  2D  and  3D  collection  scenarios  over  two 
selected  models. 

Throughout  this  study  a  control  model  was  input  to  determine  the  limitations  of  any  given 
technique.  “Model  3t”  was  created  as  a  fully-populated  noise-free  version  of  “Model  3” 
at  5cm  intervals.  For  this  model,  all  grid  and  lattice  nodes  were  filled  with  synthetic  data 
without  noise.  This  included  1681  data  points  (41x41  grid  points)  for  the  2D  data  set  and 
18,491  data  points  (41x41x1 1  lattice  points)  for  the  3D  data  set.  This  was  intended  to 
provide  perfect  input  data  to  the  algorithms  with  sufficient  resolution  to  obviate  the  need 
for  gridding,  filtering  or  other  processing.  These  results  are  shown  in  red  in  all  tables  and 
are  not  included  in  the  average  error  calculations.  A  similar  “Model  9t”  was  created 
specifically  for  the  3D  ellipsoid  inversion  and  is  covered  in  detail  in  the  appropriate 
section. 
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Figure  1:  Sample  of  3D  processing  stream  from  (a)  “cloud”  of  collected  data  points  to  (b) 
interpolated  grid  to  (c)  filtered  grid  to  (d)  analytic  signal.  View  is  from  beneath  in  the  SE  in  order  to 
illustrate  the  greater  dynamic  range  data  at  the  lowest  altitudes.  Data  shown  represent  Model  3. 
Scale  of  data  box  is  2m  x  2m  x  0.5m. 
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Figure  2:  Visualization  of  data  collection  scenarios  and  selected  models.  The  data  points  (red)  are 
common  to  analysis  of  all  models,  (a)  2D  data  collection  based  on  0.5m  line  spacing  with  random 
noise  at  slightly  varying  line  heights,  presented  over  spherical  targets  (60mm  blue,  20mm  green) 
from  Model  3.  Scale  of  box  is  approximately  2m  x  2m  x  0.6m.  (b)  3D  data  collection  based  on 
random  points  centered  over  the  target  with  noise  added,  presented  over  ellipsoid  targets  from 
Model  9.  Scale  of  box  approximately  1.5m  x  1.5m  x  0.7m. 
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Gridding,  Filtering  and  Processing 

The  2D  gridding  was  achieved  using  the  Geosoft  minimum  curvature  routine. 

Experience  has  shown  this  to  be  an  acceptable  default  for  potential  field  data. 
Representations  of  the  2D  gridded  data  are  shown  in  Appendix  B-l. 

No  such  standard  existed  for  3D  gridding.  Prior  to  choosing  a  3D  gridding  algorithm  and 
attendant  parameters,  an  investigation  was  conducted  into  the  relative  outputs  of  the 
techniques  available.  Matlab  and  Rockware  both  offered  3D  gridding  routines,  but 
Matlab  offered  a  far  more  limited  range  of  techniques.  It  was  therefore  decided  to 
concentrate  on  the  Rockware  software. 

Two  closely  spaced  objects  were  modeled  at  0.  lm  increments  in  (x,  y,  z)  and  subjected  to 
six  3D  gridding  algorithms  in  Rockworks2002.  The  six  gridding  routines  examined  were 
labeled:  closest  point,  isotropic  distance,  anisotropic  distance,  weighted  distance, 
directional  weighted  and  vertical  biased.  Each  of  these  techniques  is  described  in 
Appendix  B-5,  which  is  taken  from  the  software  help  manual.  The  directionally 
weighted  and  vertically  biased  routines  were  designed  to  emphasize  horizontal 
stratigraphic  units  only.  In  order  to  accommodate  this  orientation,  the  input  data  were 
rotated  before  gridding  so  that  up  was  North. 

Solid  models  (lattices)  were  created  with  each  algorithm  to  test  how  well  each 
represented  the  true  data.  The  quality  of  the  gridding  was  determined  as  the  standard 
deviation  of  the  difference  between  the  solid  model  output  and  the  original  modeled  data. 
The  quality  of  the  fit  between  the  lattice  results  and  the  true  model  deteriorated  with  the 
amount  of  interpolation  required  by  successive  decimation  of  the  input  data,  especially  in 
the  high  gradient  area  near  the  center  of  the  target.  Results  are  provided  in  Table  4. 

All  gridding  routines  used  a  “declustering”  technique  to  reduce  memory  requirements 
from  oversampling.  This  process  averaged  all  data  points  within  a  single  lattice  cell 
before  interpolating.  This  served  to  reduce  noise  in  closely  spaced  readings  and  to  reduce 
the  volume  of  input  to  the  gridding  routine.  None  of  these  tests  required  declustering  at 
this  stage  since  all  input  were  undersampled  rather  than  oversampled.  It  was  later  learned 
that  this  option  altered  the  position  of  the  data  by  rounding  the  position  to  the  nearest 
node  and  averaging.  This  simplified  the  interpolation  process,  but  reduced  the  intended 
accuracy  of  the  positioning  to  that  of  the  lattice  interval.  This  option  was  turned  off  for 
subsequent  operations. 

The  first  test  of  gridding  routines  used  a  fully  populated  input  data  set.  For  each  of  the 
nodes  in  the  final  lattice,  an  input  of  modeled  data  was  provided.  Under  these  conditions, 
no  interpolation  was  required  and  all  routines  faithfully  returned  identical  values  to  the 
true  modeled  data. 

The  second  test  required  interpolation  between  regularly  spaced  inputs.  Input  data  were 
provided  at  0.2m  intervals  in  (x,  y,  z)  and  were  interpolated  to  0.  lm  intervals.  This 
represented  a  data  set  desampled  to  l/8th  the  original  density.  The  best  gridding 
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algorithm  for  this  case  was  the  anisotropic,  followed  by  the  weighted  distance  and 
vertical  bias  techniques,  which  are  variants  of  the  anisotropic. 


The  third  test  also  used  l/8th  input  data  density,  but  positions  were  distributed  randomly 
throughout  the  3D  space.  The  best  algorithm  was  again  the  anisotropic  technique, 
although  the  fit  was  somewhat  poorer  than  with  the  regular  distribution.  Only  the  closest 
point  techniques  showed  a  better  fit  than  the  previous  regular  distribution. 

The  fourth  test  used  l/8th  input  data  density,  but  input  points  were  drawn  preferentially 
from  the  center  of  the  3D  space.  Half  of  the  points  were  located  randomly  within  lm  of 
the  center  point  of  the  space.  This  represents  desampling  to  1  point  out  of  every  2  within 
this  hemisphere.  The  remaining  points  were  randomly  distributed  throughout  the 
remaining  space,  representing  desampling  to  1  out  of  every  14  nodes.  All  techniques 
showed  improvement  over  the  random  distribution  due  to  the  relative  increase  in 
sampling  in  the  highest  gradient  area.  The  anisotropic  technique  again  produced  the  best 
fit. 


Table  4:  Standard  deviation  of  error  between  lattice  output  and  true  data 


Gridding  method 

closest 

point 

iso¬ 

tropic 

aniso¬ 

tropic 

dist 

wt’d 

dir 

wt’d 

vert 

wt’d 

Fully  populated 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

Desampled  points 

0.51 

0.68 

0.20 

0.25 

0.38 

0.39 

Random  distribution 

0.45 

0.69 

0.32 

0.45 

0.50 

0.45 

Concentrated  distribution 

0.41 

0.62 

0.28 

0.30 

0.31 

0.40 

Having  established  the  optimal  3D  gridding  routine  for  use  on  this  project,  filters  and 
processing  techniques  were  developed.  The  2D  filtering  and  processing  applications 
were  conducted  with  Geosoft  routines.  New  filtering  routines  had  to  be  developed  for  the 
Rockware  3D  models.  The  only  filter  available  was  an  average  filter.  This  filter  is  a 
moving  filter  which  takes  a  simple  average  of  all  points  touching  the  center  point  of  the 
moving  window  (square  roll-off).  Coefficients  for  a  3D  Hanning  filter  (cosine  roll-off) 
were  derived  and  programmed.  The  coefficients  for  the  3  point  (3x3x3  pts)  3D  Hanning 
filter  are  provided  in  Table  5. 


Table  5:  Coefficients  of  3-point  3D  Hanning  filter 


Lower  layer 

middle  layer 

upper  layer 

In  addition,  programs  were  written  to  calculate  the  analytic  signal  (AS)  from  the  3D  total 
field  data.  This  was  relatively  simple  in  that  the  vertical  gradient  could  be  calculated 
directly  from  adjacent  nodes  in  the  interpolated  lattice  just  as  the  horizontal  gradients  are 
in  2D  methods.  This  removed  the  necessity  of  using  FFT  routines.  The  Hanning  filter 
and  native  Rockware  utilities  were  all  compatible  with  the  new  AS  file  output. 
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Visualization  and  Gradient  Strings 

A  variety  of  visualization  techniques  were  available  through  the  Rockware  software  for 
the  total  field  and  analytic  signal  data.  These  included  slices  in  various  directions 
(horizontally  to  create  grids,  vertically  to  create  cross  sections,  multiple  cross  sections  to 
create  fences),  iso-surfaces  (3D  analog  of  a  single  contour  line)  and  transparencies. 

Figure  3  illustrates  some  of  these  concepts.  Color  was  used  in  all  cases  to  represent  the 
amplitude  of  the  geophysical  data.  Although  all  of  these  techniques  are  easily  understood 
and  applied  to  a  wide  variety  of  three-dimensional  data  sets,  they  require  considerable 
user  control  of  placement,  color  and  thresholds  in  order  to  adequately  visualize  the  data. 
For  the  purposes  of  smooth  potential  fields  (or  electromagnetic  fields),  a  simpler 
technique  has  been  devised  which  we  call  Gradient  Strings.  These  are  essentially  lines 
through  the  lattice  tracing  maximum  rate  of  change  in  the  field.  For  magnetic  data,  one 
can  imagine  layers  of  stacked  analytic  signal  grids  with  a  line  connecting  the  peaks  from 
one  layer  to  the  next  (Figure  4). 

The  aim  of  this  technique  is  to  reduce  the  data  set  to  its  most  important  qualities  and 
allow  simple  visualization  with  as  few  control  parameters  as  possible.  The  underlying 
assumption  is  that  maximum  gradient  lines  for  legitimate  sources  will  be  continuous  and 
will  increase  in  magnitude  with  proximity  to  the  ground.  Noise,  which  might  otherwise 
be  picked  as  a  possible  target,  will  lack  this  continuity  through  3D  space.  Similarly, 
multiple  targets  that  may  appear  as  a  single  target  on  the  original  survey  data  can  be 
resolved  as  distinct  gradient  strings.  Processing  options  based  on  Gradient  String 
information  can  include  spatial  windows,  cut-off  thresholds,  and  filter  strengths.  Strings 
can  be  down-selected  based  on  amplitude,  length  and  continuity.  If  the  line  branches 
with  depth,  the  operator  is  faced  with  the  possibility  of  multiple  targets  or  a  single  large 
target  in  the  very  near  surface. 

For  a  string  to  be  valid,  therefore,  it  must  (a)  exceed  an  intensity  threshold  (b)  be 
continuous  over  several  layers  and  (c)  increase  in  intensity  with  depth.  The  result  is  a 
line  through  the  3D  data  space  that  points  directly  to  the  target.  To  further  aid  in 
interpretation,  the  line  is  color-coded  to  gradient  intensity. 

These  gradient  strings  and  the  numerical  inversion  provide  mutual  feedback  into  the 
decision  making  process.  The  visualization  and  manipulation  of  gradient  strings  will 
allow  the  operator  to  refine  processing  parameters  and  to  deselect  targets  from  noise  or  to 
pin-point  multiple  targets  within  close  proximity  to  one  another.  In  return,  the  inversion 
will  provide  quantitative  results  in  the  form  of  fit  coefficients. 
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Figure  3:  These  graphics  illustrate  some  of  the  problems  and  complexities  associated  with 
visualization  of  3D  data,  (a)  A  solid  model  showing  the  outer  surfaces,  (b)  Slices  through  the  solid  in 
three  orthogonal  directions,  (c)  Three  iso-surfaces  at  49997nT,  50000nT  and  50017nT.  (d) 
Transparencies  of  the  iso-surfaces  combined  with  a  horizontal  slice.  Total  field  units  in  nT.  An 
example  of  a  visualization  limitation  would  be  a  poor  choice  of  iso-surface  in  (c)  which  either  missed 
or  grouped  the  anomalies. 
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Figure  4:  Illustration  of  Gradient  String  creation  from  3D  analytic  signal  of  Model  3t  (noise-free 
version  of  Model  3).  Horizontal  scale  covers  2m  x  2m,  vertical  scale  has  been  stretched  and  features 
layers  at  0.1m,  0.30m,  and  0.55m.  Analytic  signal  units  in  nT/m. 
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Detection 

For  the  purposes  of  this  study,  detection  was  defined  as  the  ability  to  reduce  the  data  to  a 
single  point  in  (x,  y)  for  each  synthetic  target.  This  basic  level  of  detection  was 
conducted  using  the  Analytic  Signal  peaks  for  the  2D  data,  and  the  Gradient  Strings  for 
the  3D  data.  Both  of  these  processes  used  the  Analytic  Signal  data  derived  from  the  total 
field  and  output  a  point  location  in  two  dimensions  that  approximated  the  center  of  the 
target.  A  10cm  grid  and  lattice  interval  was  used  throughout.  Although  the  accuracy  of 
the  Analytic  Signal  is  limited  to  the  grid  cell  size,  it  is  by  far  the  simplest  and  most 
reliable  method  of  anomaly  detection.  Subsequent  depth  and  size  calculations  (Euler, 
dipole  inversion,  etc)  may  use  total  field  data  but  the  basic  detection  is  best  achieved 
through  Analytic  Signal. 

Detection  calculations  based  on  unfiltered  data  produced  a  number  of  false  positive 
responses  with  amplitudes  comparable  to  the  primary  peaks.  A  3 -point  Hanning  filter 
was  therefore  applied  to  the  total  field  data  prior  to  calculation  of  the  Analytic  Signal  in 
order  to  smooth  the  interpolation  errors.  This  eliminated  all  false  positive  responses  from 
the  detection  list.  Smoothing  also  produced  one  false  negative  response,  missing  the 
60mm  target  in  Model  9  (Table  6).  In  this  case,  the  total  field  was  flat  enough  that  there 
was  no  Analytic  Signal  peak  above  0.  InT/m.  No  other  model  had  a  peak  less  than 
22nT/m. 

Errors  were  measured  as  the  horizontal  offset  from  the  detection  point  to  the  true  target 
location.  The  2D  data  showed  errors  that  were  slightly  higher  than  the  3D  data,  and  are 
largely  the  result  of  the  0.5m  line  spacing.  All  models,  except  for  the  fully-populated 
noise-free  version  of  Model  3  (designated  Model  3t),  picked  targets  directly  under  the 
nominal  survey  line.  This  would  indicate  that  the  introduction  of  noise  to  the  positioning 
and  magnetic  data  surpassed  the  capability  of  the  data  to  interpolate  target  locations 
between  lines,  even  at  this  relatively  tight  spacing.  The  higher  lateral  spatial  resolution  in 
the  3D  data  was  sufficient  to  overcome  this  problem. 


Table  6:  2D  Analytic  Signal  target  locations.  Average  horizontal  error  0.18m  with  one  false 
negative. 


60mm 

xy-error 

20mm 

xy-error 

True 

-0.1  -0.2 

0.1  0.3 

Model  1 

-0.0  -0.3 

0.14 

Not  modeled 

— 

Model  2 

-0.0  -0.3 

0.14 

Not  modeled 

— 

Model  3 

-0.0  -0.3 

0.14 

0.0  0.2 

0.14 

Model  3t 

-0.1  -0.3 

0.10 

0.1  0.2 

0.10 

Model  4 

-0.0  -0.3 

0.14 

Not  modeled 

— 

Model  5 

-0.0  -0.5 

0.32 

Not  modeled 

— 

Model  6 

-0.0  -0.3 

0.14 

Not  modeled 

— 

Model  7 

-0.0  -0.5 

0.32 

0.0  0.5 

0.22 

Model  8 

-0.0  -0.3 

0.14 

0.0  0.2 

0.14 

Model  9 

No  detection 

— 

0.0  0.2 

0.14 
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As  described  by  Salem  et  al.  (2002),  there  is  a  consistent  shift  of  the  AS  peak  to  the 
magnetic  South  by  an  amount  dependent  upon  the  target  depth  and  the  field  inclination. 
Without  knowing  the  true  depth,  this  cannot  be  calculated  to  adjust  the  AS  peak 
locations.  The  average  shift  observed,  however,  is  consistent  with  the  typical  range  of 
UXO  burial  depths.  As  will  be  discussed  later,  the  use  of  the  depth  estimates  from  other 
techniques  can  be  used  to  refine  the  horizontal  position  estimate  of  the  Analytic  Signal. 

From  the  filtered  3D  data,  Gradient  Strings  were  calculated  through  the  solid  lattice  of 
the  Analytic  Signal  data.  The  most  consistent  results  were  obtained  using  a  10cm  lattice 
size.  Smaller  intervals  displayed  sufficient  noise  and  non-dipole  effects  to  disconnect  the 
Gradient  Strings.  The  larger  grid  smoothed  the  noise  and  allowed  the  gradient  to  be 
calculated  over  a  longer  baseline.  The  10cm  baseline  was  generally  longer  than  the 
wavelength  of  the  non-dipole  effects,  thus  effectively  masking  them.  Strings  less  than  3 
points  in  length,  and  having  a  maximum  value  <35nT/m  were  deleted  as  noise.  In  all 
cases,  the  remaining  strings  were  continuously  increasing  in  amplitude  with  proximity  to 
the  ground.  Graphical  results  are  shown  in  Appendix  B-2. 

Most  of  the  models  with  the  deeply  buried  60mm  target  (2, 3, 6, 7, 8)  had  difficulty 
maintaining  a  continuous  string  throughout  the  entire  data  volume.  Usually,  strong 
gradient  noise  at  an  intermediate  height  split  the  strand  into  two  discontinuous  parts  and 
had  to  be  reconnected  manually.  Models  4,  5  and  9  had  the  same  target  depth  but  were 
able  to  maintain  a  continuous  string.  It  may  therefore  be  concluded  that  the  Gradient 
String  method  has  some  limits  regarding  its  ability  to  maintain  a  continuous  string. 

These  limits  appear  to  be  reached  only  at  the  maximum  expected  burial  depth  with  the 
moderately  high  noise  levels  provided  and  can  be  solved  with  manual  inspection.  Further 
experimentation  can  be  expected  to  remedy  this  minor  flaw  in  automation. 

Target  location  with  this  method  was  within  one  lattice  point  in  every  case  but  one.  A 
smaller  lattice  size  improved  resolution  but  decreased  string  continuity.  Final  results  for 
each  model  are  presented  in  Table  7.  Numbers  in  brackets  are  averages  over  the  length 
of  the  string. 


Table  7:  Gradient  String  target  locations.  Average  horizontal  error  0.09m. 


60mm 

xy-error 

20mm 

xy-error 

True 

-0.1  -0.2 

0.1  0.3 

Model  1 

-0.2  -0.2  (-0.12-0.28) 

0.1  (0.08) 

Not  modeled 

— 

Model  2 

-0.1  -0.1  (-0.13  -0.28) 

0.1  (0.09) 

Not  modeled 

— 

Model  3 

-0.2  -0.2  (-0.20  -0.33) 

0.1 (0.17) 

0.0  0.3  (0.08  0.27) 

0.1  (0.04) 

Model  3t 

-0.1  -0.3  (-0.10-0.32) 

0.1 (0.12) 

0.1  0.3  (0.10  0.26) 

0.0  (0.04) 

Model  4 

-0.2  -0.2  (-0.15  -0.32) 

0.1  (0.13) 

Not  modeled 

— 

Model  5 

-0.2  -0.2  (-0.18-0.27) 

0.1  (0.11) 

Not  modeled 

— 

Model  6 

-0.1  -0.3  (-0.12-0.22) 

0.1  (0.02) 

Not  modeled 

— 

Model  7 

-0.2  -0.2  (-0.18  -0.28) 

0.1  (0.12) 

0.0  0.3  (0.08  0.25) 

0.1  (0.06) 

Model  8 

-0.1  -0.3  (-0.14-0.30) 

0.1  (0.11) 

0.1  0.4(0.10  0.32) 

0.1  (0.02) 

Model  9 

-0.0  -0.4  (-0.12-0.30) 

0.4(0.10) 

0.1  0.3  (0.08  0.23) 

0.0  (0.08) 
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In  conclusion,  the  gradient  string  method  successfully  identified  the  number  of  targets  in 
each  model  with  an  automated  approach,  and  located  them  within  one  grid  cell  of  the  true 
location.  A  10cm  lattice  was  used  to  reduce  the  effects  of  strong  gradients  in  the  analytic 
signal  at  the  lowest  heights,  which  affected  the  final  target  location  accuracy.  Some 
models  required  visual  inspection  to  connect  a  pair  of  disjointed  gradient  strings.  It  is 
worth  noting  that  the  inaccuracy  in  the  target  location  is  partially  the  result  of  the 
theoretical  limitations  of  the  analytic  signal,  as  demonstrated  by  the  noise-free  data.  This 
can  be  reduced  by  feedback  using  depth  estimates  from  other  analysis  techniques. 
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Depth  Estimation 

Calculations  for  target  depth  were  conducted  using  the  ratio  of  analytic  signals  (after 
Salem  et  ah,  2002).  This  technique  uses  the  equation 

Z=4*ASi/AS2 

where  ASi  is  the  analytic  signal  of  the  total  field  and  AS2  is  the  analytic  signal  of  the 
vertical  gradient  (a  second  order  gradient,  units  of  nT/m2). 

For  the  2D  data,  these  parameters  were  calculated  using  the  Geosoft  FFT  algorithms  from 
the  gridded  data.  For  the  3D  data,  they  were  calculated  from  the  total  field  lattice  data. 
The  point  of  calculation  is  the  peak  of  the  ASi.  For  the  3D  data,  this  was  taken  at  the 
average  (x,  y)  location  of  the  ASi  peak,  but  at  a  height  of  0.3m.  This  accounted  for  the 
boundary  errors  inherent  in  the  calculation.  For  example,  if  the  data  at  0.  lm  were  used, 
the  vertical  gradient  calculations  would  be  based  on  data  outside  the  bounds  of  the 
collected  data  space.  It  also  made  the  2D  and  3D  results  more  comparable,  since  they 
were  calculated  at  0.3m  and  0.35m  respectively. 

It  was  also  found  that  the  unfiltered  total  field  data  produced  significantly  better  depth 
estimates  than  the  filtered  data.  The  application  of  a  filter  reduced  and  broadened  the 
peak  amplitudes  thereby  increasing  the  apparent  depth  of  the  target  by  approximately  two 
times.  Final  depths  and  errors  are  presented  in  Tables  8  &  9.  In  most  cases,  the  vertical 
accuracy  is  the  same  or  better  than  the  horizontal  accuracy.  The  generally  poorer 
performance  of  the  20mm  depth  estimates  is  attributed  to  the  lower  signal-noise  ratio  at 
the  0.3m  calculation  height.  The  use  of  a  smaller  lattice  interval  may  have  improved 
these  results  by  reducing  the  effective  calculation  height  thereby  increasing  the  signal- 
noise  ratio  at  the  point  of  calculation.  Instead,  it  was  decided  to  maintain  the  1 0cm 
interval  from  the  Gradient  String  calculation  for  the  sake  of  continuity  and  comparison  to 
the  2D  data. 

The  2D  depth  estimates  showed  nearly  five  times  greater  error  than  the  3D  data, 
averaging  0.29m  depth  error  for  the  60mm  compared  to  0.06m  for  the  3D  data.  In 
general,  the  depth  of  the  60mm  is  underestimated  by  one-half  while  the  20mm  is 
overestimated.  It  is  interesting  to  note  that  the  noise-free  data  accurately  returned  the 
depth  of  the  60mm,  but  not  the  20mm  target. 

The  3D  data  showed  exceptionally  good  depth  estimates  for  the  60mm,  but  did  less  well 
on  the  smaller  20mm  target.  Having  calculated  the  depth,  it  would  now  be  possible  to 
refine  the  horizontal  positions  from  the  Gradient  Strings  and  further  improve  the  accuracy 
to  approximately  5cm  radially  from  the  target  center.  This  is  an  iterative  feedback 
approach  that  should  be  explored  in  more  detail  in  subsequent  studies. 
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Table  8:  2D  analytic  signal  depth  estimates.  Average  depth  error  0.21m. 


60mm 

z-error 

20mm 

z-error 

True 

0.60* 

0.05 

Model  1 

0.10 

0.10 

Not  modeled 

— 

Model  2 

0.22 

0.38 

Not  modeled 

— 

Model  3 

0.29 

0.31 

0.04 

0.01 

Model  3t 

0.63 

0.03 

0.13 

0.08 

Model  4 

0.03 

0.17 

Not  modeled 

— 

Model  5 

0.22 

0.38 

Not  modeled 

— 

Model  6 

0.30 

0.30 

Not  modeled 

— 

Model  7 

0.23 

0.37 

0.02 

0.03 

Model  8 

0.37 

0.23 

0.10 

0.05 

Model  9 

0.2 1 1 

0.39 

0.08 

0.03 

*depth  of  models  1  and  4  is  0.20m 

twhereas  no  peak  was  detected  for  this  target  in  the  fdtered  data,  the  depths  were 
calculated  from  unfdtered  data  which  did  show  a  peak  at  location  (-0.0,  -0.5) 


Table  9:  3D  analytic  signal  depth  estimates.  Average  depth  error  0.10m. 


60mm 

z-error 

20mm 

z-error 

True 

0.60* 

0.05 

Model  1 

0.07 

0.13 

Not  modeled 

— 

Model  2 

0.56 

0.04 

Not  modeled 

— 

Model  3 

0.59 

0.01 

0.25 

0.20 

Model  3t 

0.57 

0.03 

0.00 

0.05 

Model  4 

0.23 

0.03 

Not  modeled 

— 

Model  5 

0.56 

0.04 

Not  modeled 

— 

Model  6 

0.58 

0.02 

Not  modeled 

— 

Model  7 

0.58 

0.02 

0.01 

0.04 

Model  8 

0.48 

0.12 

0.27 

0.23 

Model  9 

0.47 

0.13 

0.29 

0.24 

*depth  of  models  1  and  4  is  0.20m 
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Euler  Deconvolution 

The  Euler  equation 


{x~xo)~~Jr  (y  ~  Vq)~~  +  (z  ~  zo)4~  =  N(F  —  T) 
ox  oy  Oz 


can  be  written  as 


dT  dT  dT  dT  dT  dT  „ 

— +-z-^0  +^zo  =^x+^y+^z-N(F~T) 

ox  oy  oz  ox  oy  oz 


or 


x0 

dT  dT  dT 

* 

dx  dy  dz 

To 

_v 

-M 


where  (x,  y,  z)  is  the  sample  location,  (xo,  yo,  zo)  is  the  target  location,  T  is  the  measured 
total  field  around  the  anomaly,  F  is  the  background  magnetic  field,  and  N  is  an  index 
number  associated  with  the  dimensionality  of  the  target. 

A  window  of  points  around  the  sample  location  produces  a  set  of  equations  of  this  fonn 
which  can  be  solved  as  a  series  of  linear  equations  in  the  form  AX=B.  The  size  of  the 
window  is  generally  determined  by  the  spatial  extent  of  the  anomaly  in  order  to  improve 
the  signal  content  of  the  least  squares  equations.  This  window  should  be  large  enough  to 
encompass  a  reasonable  portion  of  the  anomaly,  but  small  enough  that  it  does  not  include 
signals  from  overlapping  anomalies.  For  these  tests,  a  simple  3x3  or  3x3x3  window  was 
used  to  produce  9  or  27  equations  for  the  2D  and  3D  respectively.  The  index  parameter 
N  was  set  to  2.8  to  achieve  the  best  results  in  both  2D  and  3D  processing.  As  with  the 
analytic  signal  depth  calculations,  better  results  were  obtained  using  the  unfiltered  rather 
that  the  filtered  total  field  data. 

As  with  the  analytic  signal  depth  calculation,  the  Euler  calculation  was  run  at  the  average 
(x,  y)  location  specified  by  the  gradient  calculation  and  at  a  height  of  0.3m.  This  allowed 
a  buffer  zone  for  the  vertical  gradient  calculation  and  conveniently  corresponded  to  the 
height  of  the  2D  calculation  for  closer  comparison. 

Euler  deconvolution  of  the  2D  data  was  conducted  using  the  Geosoft  grid  routines,  since 
the  intention  was  to  compare  the  3D  method  against  current  industry  standards.  These 
are  more  sophisticated  in  their  data  conditioning  than  was  possible  for  the  3D  data  in  this 
project.  This  condition  notwithstanding,  the  3D  technique  again  produced  average  depth 
errors  approximately  five  times  smaller  than  the  2D  data  (0.03m  and  0.14m  respectively 
for  the  60mm  target).  This  is  roughly  half  that  of  the  AS  ratio  technique  described  above. 
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Table  10:  2D  Euler  depth  estimates.  Average  depth  error  0.24m. 


60mm 

z-error 

20mm 

z-error 

True 

0.60* 

0.05 

Model  1 

0.23  J 

0.03 

Not  modeled 

— 

Model  2 

0.46 

0.14 

Not  modeled 

— 

Model  3 

0.45 

0.15 

0.56 

0.51 

Model  3t 

0.57 

0.03 

0.24 

0.19 

Model  4 

0.23J 

0.03 

Not  modeled 

— 

Model  5 

0.58 

0.02 

Not  modeled 

— 

Model  6 

0.26 

0.34 

Not  modeled 

— 

Model  7 

0.15 

0.45 

0.43 

0.38 

Model  8 

0.48 

0.12 

0.53 

0.48 

Model  9 

0.49| 

0.11 

0.47 

0.42 

*depth  of  models  1  and  4  is  0.20m 
^required  manual  adjustment  of  window  size 

twhereas  no  peak  was  detected  for  this  target  in  the  filtered  data,  the  depths  were 
calculated  from  unfiltered  data  which  did  show  a  peak  at  location  (-0.0,  -0.5) 


Table  11:  3D  Euler  depth  estimates.  Average  depth  error  0.06m. 


60mm 

z-error 

20mm 

z-error 

True 

0.60* 

0.05 

Model  1 

0.19 

0.01 

Not  modeled 

— 

Model  2 

0.61 

0.01 

Not  modeled 

— 

Model  3 

0.64 

0.04 

0.19 

0.14 

Model  3t 

0.61 

0.01 

0.03 

0.02 

Model  4 

0.24 

0.04 

Not  modeled 

— 

Model  5 

0.56 

0.04 

Not  modeled 

— 

Model  6 

0.57 

0.03 

Not  modeled 

— 

Model  7 

0.58 

0.02 

0.07 

0.04 

Model  8 

0.63 

0.03 

0.28 

0.23 

Model  9 

0.62 

0.02 

0.20 

0.15 

*depth  of  models  1  and  4  is  0.20m 
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Dipole  Inversion 

The  first  inversion  for  parameters  other  than  depth  was  conducted  using  the  MTADS- 
DAS  dipole  inversion  code.  This  code  operates  on  a  set  of  (x,  y,  z,  mag)  data  that  may  or 
may  not  be  broken  into  lines.  As  such  it  is  well  suited  to  this  project.  It  uses  its  own 
gridding  algorithm  to  create  a  presentation  data  set  that  is  used  to  manually  isolate 
anomalies  for  inversion,  but  it  operates  on  the  original  point  data  set.  It  is  therefore 
independent  of  the  2D  or  3D  gridding  process  and  makes  full  use  of  3D  positioning  data, 
but  does  not  benefit  from  the  advantages  of  down-line,  cross-line,  or  vertical  smoothing 
through  interpolation  or  filtering. 

This  code  operates  on  one  anomaly  at  a  time.  Where  signatures  overlap,  data  must  be 
manually  divided  between  anomalies  based  on  position.  For  this  project,  where  the 
detection  algorithms  identified  multiple  targets  (Models  3,  7,  8,  9)  the  dividing  line  was 
drawn  from  (x,  y)  coordinates  (-1,  l)-(0,  0)-(l,  0).  Figure  5  illustrates  the  demarcation 
line  used  for  segmentation.  This  effectively  bisected  the  total  field  peaks,  but  the  troughs 
were  too  closely  overlapped  to  separate  and  were  therefore  attributed  to  the  more 
northerly  20mm  target.  This  segmentation  process  made  it  difficult  for  the  program  to 
invert  the  data  in  several  instances.  For  this  reason,  where  multiple  anomalies  were 
detected,  they  were  also  inverted  as  a  single  “combined”  anomaly  for  comparison.  This 
is  comparable  to  the  UXAnalyze  inversion  results  discussed  below. 

The  most  obvious  difference  between  the  2D  and  3D  data  was  the  failure  of  the  code  to 
invert  the  2D  data  for  any  of  the  20mm  targets  and  half  of  the  60mm  targets.  The  DAS 
code  inverted  the  3D  data  for  all  targets  except  for  Model  6.  There  were  also  some 
irregularities  in  the  fit  coefficients  that  have  not  been  explained.  For  example,  the  noise- 
free  data  in  the  2D  Model  3t  were  inverted  for  the  combined  anomalies  with  a  better  fit 
than  the  60mm  alone.  The  2D  Model  7  was  able  to  invert  the  single  60mm,  but  did  not 
have  sufficient  signal  to  invert  the  combined  anomalies.  The  3D  Model  1  was  inverted 
for  the  60mm  for  position  with  reasonable  accuracy,  but  output  a  fit  coefficient  of  0. 

The  poor  perfonnance  is  attributed  primarily  to  data  segmentation  problems,  but  there 
were  also  difficulties  with  the  IDL-based  version  of  the  MTADS-DAS  dipole  inversion 
code.  It  is  more  restrictive  than  the  UXAnalyze  code  currently  under  development  by 
AETC  and  Geosoft  which  uses  the  same  inversion  engine.  It  was  specifically  designed  to 
process  and  analyze  MTADS  survey  data,  and  as  a  result  many  parameters  are  hardwired 
or  defaulted  to  MTADS  specifications.  This  made  it  extremely  difficult  to  implement  for 
this  3D  SEED  Project.  AETC  generously  contributed  the  time  and  effort  to  invert  these 
Models  using  their  UXAnalyze  code  while  still  under  development.  The  results  are 
shown  in  Appendix  B-3.  This  resolved  more  of  the  2D  data  sets  than  the  IDL-based  code 
used  here,  but  where  both  codes  reached  a  solution  the  accuracy  of  the  target  positioning 
was  comparable.  Also,  if  the  anomalies  are  treated  as  “combined”  targets,  the  number  of 
inversion  failures  is  the  same  for  both  methods. 
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Figure  5:  Screen  capture  showing  demarcation  line  used  to  segment  data  into  separate  anomalies 
where  multiple  anomalies  were  detected. 
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Table  12:  2D  data  fit  coefficients. 


Fit  (60mm) 

Fit  (20mm) 

Fit  (comb) 

Model  1 

0.976 

Not  modeled 

— 

Model  2 

0.973 

Not  modeled 

— 

Model  3 

Not  strong  enough 

Not  strong  enough 

0.940 

Model  3t 

0.890 

Not  strong  enough 

0.908 

Model  4 

0.957 

Not  modeled 

— 

Model  5 

Not  strong  enough 

Not  modeled 

— 

Model  6 

0.982 

Not  modeled 

— 

Model  7 

0.949 

Not  strong  enough 

Not  strong  enough 

Model  8 

Not  strong  enough 

Not  strong  enough 

0.951 

Model  9 

Not  strong  enough 

Not  strong  enough 

0.846 

Table  13:  2D  dipole  inversion  results.  Average  location  error  0.12m  with  nine  failures. 


60mm 

Error 

20mm 

Error 

X  Y  Z  size 

X  Y  Z  size 

True 

-0.10  -0.20  0.60*  0.10 

0.10  0.30  0.05  0.03 

Model  1 

-0.18  -0.18  0.16  0.05 

0.09 

Not  modeled 

Model  2 

-0.08  -0.13  0.49  0.04 

0.13 

Not  modeled 

Model  3 

— 

— 

— 

— 

Model  3t 

-0.20  -0.13  0.31  0.03 

0.32 

— 

— 

Model  4 

-0.11  -0.18  0.20  0.05 

0.02 

Not  modeled 

Model  5 

— 

— 

Not  modeled 

Model  6 

-0.10  -0.11  0.49  0.06 

0.14 

Not  modeled 

Model  7 

-0.22  -0.20  0.45  0.04 

0.19 

— 

— 

Model  8 

— 

— 

— 

— 

Model  9 

— 

— 

— 

— 

*depth  of  models  1  and  4  is  0.20m 
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Table  14:  3D  data  fit  coefficient. 


Fit  (60mm) 

Fit  (20mm) 

Fit  (comb) 

Model  1 

0.000 

Not  modeled 

— 

Model  2 

0.791 

Not  modeled 

— 

Model  3 

0.948 

0.693 

0.299 

Model  3t 

0.913 

0.004 

0.001 

Model  4 

0.823 

Not  modeled 

— 

Model  5 

0.494 

Not  modeled 

— 

Model  6 

Too  much  change 

Not  modeled 

— 

Model  7 

0.421 

0.658 

X  optimize  prob 

Model  8 

0.922 

0.635 

0.603 

Model  9 

0.773 

0.166 

0.513 

Table  15:  3D  dipole  inversion  results.  Average  location  error  0.28m  with  one  failure. 


60mm 

Error 

20mm 

Error 

X 

Y 

Z 

size 

X 

Y 

Z 

size 

True 

-0.10 

-0.20 

0.60* 

0.10 

0.10 

0.30 

0.05 

0.03 

Model  1 

-0.05 

-0.28 

-0.12 

0.00 

0.33 

Not  modeled 

Model  2 

-0.26 

-0.28 

0.35 

0.03 

0.31 

Not  modeled 

Model  3 

-0.28 

-0.20 

0.43 

0.04 

0.25 

0.03 

0.16 

0.04 

0.01 

0.16 

Model  3t 

-0.17 

-0.15 

0.46 

0.04 

0.16 

0.10 

0.37 

-0.10 

0.00 

0.17 

Model  4 

-0.17 

-0.24 

0.14 

0.03 

0.10 

Not  modeled 

Model  5 

-0.35 

-0.22 

0.18 

0.02 

0.49 

Not  modeled 

Model  6 

— 

— 

Not  modeled 

Model  7 

-0.25 

-0.33 

0.16 

0.01 

0.48 

0.02 

0.34 

-0.01 

0.01 

0.11 

Model  8 

+0.05 

-0.32 

0.28 

0.04 

0.37 

0.03 

0.16 

0.04 

0.01 

0.16 

Model  9 

-0.15 

-0.19 

0.19 

0.03 

0.41 

0.07 

0.23 

-0.16 

0.01 

0.22 

*depth  of  models  1  and  4  is  0.20m 


The  positions  reported  by  the  DAS  inversion  results  were  comparable  in  (x,  y)  accuracy 
to  the  analytic  signal  and  Gradient  String  results.  The  largest  contribution  to  radial 
position  error  came  from  the  depth  estimates,  which  were  underestimated  in  all  cases. 
Several  depths  were  calculated  as  above  the  ground  surface.  The  60mm  depth  estimates 
from  the  3D  data  set  showed  greater  average  error  than  the  2D  data  (0.25m  and  0.12m 
respectively).  In  general,  the  dipole  inversion  of  the  2D  models  generated  better  depth 
estimates  than  the  2D  Euler  and  analytic  signal  ratio  methods,  but  the  dipole  inversion  of 
3D  models  was  4-10  times  worse.  There  were  so  few  successful  2D  models  however  that 
this  is  a  doubtful  statistic. 
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Although  inclination,  azimuth  and  moment  information  were  also  generated  by  this 
program,  these  parameters  are  linked  to  remanence  rather  than  target  orientation.  Since 
this  was  not  a  variable  in  the  target  models,  the  expectation  is  that  these  parameters 
should  align  with  the  background  field  within  +/-  30°.  The  observed  results  were 
generally  consistent  with  this  expectation,  although  results  may  be  significantly  distorted 
by  the  segmentation  process,  which  separated  much  of  the  negative  lobe  of  the  60mm 
target  from  the  rest  of  the  signature. 

There  are  at  least  three  reasons  for  the  generally  poor  performance  of  this  technique.  The 
first  is  the  requirement  for  segmentation  of  the  data.  As  has  been  stated  above,  this  had 
the  effect  of  separating  the  negative  and  positive  peaks  in  some  of  the  strongest 
anomalies.  The  second  was  the  possible  mis-match  in  attempting  to  fit  a  point  dipole  to 
an  ellipsoidal  model.  This  was  reflected  in  the  poorer  results  of  Models  4-9.  Inversion 
for  multiple  targets,  and  for  shapes  other  than  a  dipole  is  presented  in  the  next  section.  A 
third  problem  was  that  the  inclination,  declination  and  background  field  strength  were 
hardwired  to  the  IGRF  lookup  table  at  the  latitude  and  longitude  provided.  Although  the 
inclination  and  declination  could  be  adequately  represented,  the  field  strength  of 
50,000nT  could  not  be.  This  infonnation  is  not  provided  by  the  program,  but  it  is 
estimated  that  the  field  strength  used  by  the  software  was  approximately  55,000nT.  The 
program  would  therefore  expect  much  stronger  anomalies  than  the  input  data  presented, 
and  would  consequently  underestimate  the  depths  of  the  calculated  dipoles. 
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Solid-body  Inversion 

Inversion  of  the  data  for  shaped  bodies  (spheres  and  ellipsoids)  was  conducted  using  the 
ModelVision  software  from  Encom.  This  program  was  used  to  invert  for  spheres  and 
ellipsoids  from  a  pre-set  starting  model.  The  vendor  recommends  a  more  finely  tuned 
starting  model  through  manual  manipulation  of  target  parameters  for  each  data  set,  but 
this  was  ignored  in  order  to  assure  unifonnity  for  testing  purposes.  This  software  had  the 
advantage  of  inverting  for  multiple  targets  without  the  need  for  segmentation.  The 
disadvantage  is  that  the  fit  coefficient  is  an  aggregate  over  the  entire  data  set  rather  than 
being  target  specific. 

The  2D  and  3D  data  sets  were  inverted  once  for  spheres,  and  then  again  for  ellipsoids. 

For  spheres,  starting  models  included  one  or  two  targets  (as  detennined  by  the  previous 
detection  procedure)  with  locations  of  (-0.1,  -0.3,  0.5)  and  (0.1,  0.3,  0.5)  with  a  radius  of 
0.1m  for  both.  Ellipsoid  starting  models  were  centered  on  the  same  locations,  but  had 
axes  lengths  of  (major  0.05,  semi-major  0.02,  minor  0.02)  and  orientations  of  (azimuth 
90°,  plunge  0°,  rotation  0°).  Maximum  ranges  for  each  variable  were  set  to  allow 
sufficient  change  in  the  model,  and  a  maximum  of  20  iterations  were  allowed. 

Each  spherical  inversion  included  four  free  inversion  parameters:  X,  Y,  Z,  radius.  Each 
ellipsoidal  inversion  included  nine  free  inversion  parameters:  X,  Y,  Z,  axis-A,  axis-B, 
axis-C,  azimuth,  plunge,  rotation.  Susceptibility  was  fixed  at  3  SI  for  all  models.  Models 
with  two  targets  (Models  3,  7,  8,  9)  doubled  the  number  of  free  parameters.  RMS  error 
for  multiple  target  models  is  an  aggregate  fit  over  the  entire  data  set.  It  should  be  noted 
here  that  Models  1-3  were  synthesized  from  spheres,  whereas  Models  4-9  were 
synthesized  from  ellipsoids. 

Tables  16-19  show  the  2D  inversion  results,  while  Tables  20-23  show  the  comparable  3D 
inversion  results.  Tables  16  and  20  show  the  general  fit  coefficients.  Tables  17  and  21 
show  the  inversion  results  (positions)  when  testing  for  spherical  models.  Tables  18/19 
and  22/23  show  the  inversion  results  (positions  and  dimensions)  when  testing  for 
ellipsoidal  models. 
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Table  16:  2D  RMS  error  when  fitting  a  sphere  or  ellipsoid  to  the  data. 


RMS  error  (sph) 

RMS  error  (ellip) 

Model  1 

1.25 

1.23 

Model  2 

1.29 

1.28 

Model  3 

1.25 

1.47 

Model  3t 

0.01 

0.03 

Model  4 

5.52 

5.46 

Model  5 

8.30 

8.23 

Model  6 

2.49 

2.50 

Model  7 

7.50 

7.45 

Model  8 

2.40 

2.06 

Model  9 

3.83 

3.80 

Table  17:  2D  Spherical  inversion  results  (fitting  a  spherical  model).  Average  location  error  0.14m. 


60mm 

Error 

20mm 

Error 

X 

Y 

z 

rad 

X  Y 

Z 

rad 

True 

-0.10 

-0.20 

0.60* 

0.05 

0.10  0.30 

0.05 

0.02 

Model  1 

-0.10 

-0.20 

0.20 

0.05 

0.00 

Not  modeled 

— 

Model  2 

-0.10 

-0.20 

0.60 

0.05 

0.00 

Not  modeled 

— 

Model  3 

-0.10 

-0.20 

0.60 

0.05 

0.00 

-0.07  0.30 

0.07 

0.02 

0.03 

Model  3t 

-0.10 

-0.20 

0.60 

0.05 

0.00 

0.10  0.30 

0.05 

0.02 

0.00 

Model  4 

-0.11 

-0.38 

0.19 

0.05 

0.18 

Not  modeled 

— 

Model  5 

-0.08 

-0.50 

0.56 

0.04 

0.30 

Not  modeled 

— 

Model  6 

-0.11 

-0.10 

0.61 

0.07 

0.10 

Not  modeled 

— 

Model  7 

-0.11 

-0.50 

0.57 

0.04 

0.30 

0.21  0.10 

0.00 

0.02 

0.23 

Model  8 

-0.11 

-0.16 

0.59 

0.06 

0.04 

-0.02  0.34 

0.31 

0.04 

0.28 

Model  9 

-0.23 

-0.10 

0.59 

0.06 

0.17 

-0.08  0.28 

0.18 

0.03 

0.22 

*depth  of  models  1  and  4  is  0.20m 
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Table  18:  2D  Ellipsoid  inversion  results  (positions).  Average  location  error  0.17m. 


60mm 

Error 

20mm 

Error 

X 

Y 

Z 

X  Y 

Z 

True 

-0.10 

-0.20 

0.60* 

0.10  0.30 

0.05 

Model  1 

-0.10 

-0.20 

0.20 

0.00 

Not  modeled 

— 

Model  2 

-0.10 

-0.20 

0.61 

0.01 

Not  modeled 

— 

Model  3 

-0.10 

-0.23 

0.61 

0.03 

0.01  0.27 

0.18 

0.16 

Model  3t 

-0.10 

-0.20 

0.61 

0.01 

0.10  0.30 

0.05 

0.00 

Model  4 

-0.10 

-0.38 

0.21 

0.14 

Not  modeled 

— 

Model  5 

-0.09 

-0.49 

0.55 

0.29 

Not  modeled 

— 

Model  6 

-0.11 

-0.10 

0.62 

0.10 

Not  modeled 

— 

Model  7 

-0.10 

-0.50 

0.57 

0.30 

0.21  0.10 

0.00 

0.23 

Model  8 

-0.12 

-0.17 

0.58 

0.04 

-0.03  0.32 

0.43 

0.40 

Model  9 

-0.24 

-0.10 

0.58 

0.17 

-0.09  0.29 

0.21 

0.25 

*depth  of  models  1  and  4  is  0.20m 


Table  19:  2D  Ellipsoid  inversion  results  (dimensions) 


60mm 

Error 

20mm 

Error 

A 

B 

C 

%Vol 

A  B 

C 

%Vol 

True 

0.051 

0.051 

0.051 

0.017  0.017 

0.017 

Model  1 

0.10 

0.04 

0.04 

1 

Not  modeled 

— 

Model  2 

0.15 

0.04 

0.03 

2 

Not  modeled 

— 

Model  3 

0.10 

0.05 

0.03 

3 

0.04  0.02 

0.01 

104 

Model  3t 

0.10 

0.05 

0.03 

0 

0.03  0.02 

0.02 

1 

True 

0.100 

0.030 

0.030 

0.035  0.010 

0.010 

Model  4 

0.14 

0.03 

0.02 

7 

Not  modeled 

— 

Model  5 

0.13 

0.04 

0.02 

6 

Not  modeled 

— 

Model  6 

0.17 

0.07 

0.03 

256 

Not  modeled 

— 

Model  7 

0.12 

0.05 

0.02 

1 

0.05  0.01 

0.01 

42 

Model  8 

0.21 

0.05 

0.02 

170 

0.36  0.01 

0.01 

1925 

Model  9 

0.11 

0.06 

0.03 

162 

0.13  0.01 

0.01 

613 
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The  most  notable  result  of  the  2D  inversion  was  that  the  spherical  targets  (Models  1-3) 
were  better  resolved  than  the  ellipsoids.  This  was  seen  most  clearly  in  the  RMS  errors.  It 
should  also  be  noted  that  the  z-error  was  no  longer  the  dominant  position  error. 

Horizontal  errors  in  target  location,  especially  for  the  ellipsoid  targets,  were  larger  than 
those  previously  detennined  by  the  gradient  strings.  The  location  error  presented  here  is 
the  total  offset  including  both  horizontal  and  vertical  components. 

The  high  errors  in  Models  4-9  were  attributed  to  the  differences  in  ellipsoid  modeling 
techniques  between  MAGMOD  and  ModelVision.  Codes  for  spherical  models  compare 
exceptionally  well  (Models  1-3),  but  as  has  already  been  stated,  these  codes  produce 
significantly  different  results  for  ellipsoidal  bodies.  These  effects  can  only  be  attributed 
to  the  difference  in  approach  to  modeling  between  MAGMOD  and  ModelVision.  This  is 
further  explored  in  the  noise-free  Model  9t  results  presented  below. 

Ellipsoid  inversion  results  for  target  orientation  were  calculated  but  were  not  explored. 
The  ModelVision  software  inverts  for  full  ellipsoids,  and  not  cylindrical  ellipsoids  -  it 
does  not  allow  the  variation  in  the  lesser  axes  to  be  tied  together.  As  a  result,  inversion 
targets  may  appear  as  disks  or  ovals  rather  than  round  cylinders,  which  made  orientation 
comparisons  difficult.  This,  added  to  the  large  errors  in  the  target  dimensions,  made  the 
orientation  analysis  irrelevant. 
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Table  20:  3D  RMS  error  when  fitting  a  sphere(s)  or  ellipsoid(s)  to  the  data. 


RMS  error  (sph) 

RMS  error  (ellip) 

Model  1 

1.29 

1.12 

Model  2 

1.29 

1.29 

Model  3 

0.74 

0.71 

Model  3t 

0.01 

0.02 

Model  4 

3.85 

3.94 

Model  5 

6.23 

6.26 

Model  6 

1.65 

1.75 

Model  7 

3.52 

3.51 

Model  8 

1.44 

6.55 

Model  9 

1.25 

1.27 

Model  9t 

1.09 

1.41 

Table  21:  3D  Spherical  inversion  results  (fitting  a  spherical  model).  Average  location  error  0.09m. 


60mm 

Error 

20mm 

Error 

X 

Y 

Z 

rad 

X 

Y 

Z 

rad 

True 

-0.10 

-0.20 

0.60* 

0.05 

0.10 

0.30 

0.05 

0.02 

Model  1 

-0.10 

-0.21 

0.21 

0.05 

0.01 

Not  modeled 

— 

Model  2 

-0.10 

-0.20 

0.60 

0.05 

0.00 

Not  modeled 

— 

Model  3 

-0.10 

-0.20 

0.60 

0.05 

0.01 

0.11 

0.30 

0.04 

0.02 

0.01 

Model  3t 

-0.10 

-0.20 

0.60 

0.05 

0.00 

0.10 

0.30 

0.05 

0.02 

0.00 

Model  4 

-0.07 

-0.32 

0.17 

0.04 

0.13 

Not  modeled 

— 

Model  5 

-0.07 

-0.47 

0.50 

0.04 

0.29 

Not  modeled 

— 

Model  6 

-0.10 

-0.13 

0.62 

0.07 

0.08 

Not  modeled 

— 

Model  7 

-0.08 

-0.47 

0.56 

0.05 

0.27 

0.12 

0.22 

0.03 

0.01 

0.08 

Model  8 

-0.11 

-0.14 

0.63 

0.07 

0.07 

0.10 

0.33 

0.05 

0.02 

0.03 

Model  9 

-0.18 

-0.10 

0.60 

0.06 

0.13 

0.04 

0.29 

0.06 

0.02 

0.06 

Model  9t 

-0.19 

-0.10 

0.60 

0.06 

0.13 

0.04 

0.29 

0.05 

0.02 

0.06 

*depth  of  models  1  and  4  is  0.20m 
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Table  22:  3D  Ellipsoid  inversion  results  (positions).  Average  location  error  0.18m. 


60mm 

Error 

20mm 

Error 

X 

Y 

Z 

X  Y 

Z 

True 

-0.10 

-0.20 

0.60* 

0.10  0.30 

0.05 

Model  1 

-0.10 

-0.20 

0.19 

0.01 

Not  modeled 

— 

Model  2 

-0.10 

-0.20 

0.62 

0.02 

Not  modeled 

— 

Model  3 

-0.10 

-0.19 

0.60 

0.01 

0.10  0.30 

0.05 

0.01 

Model  3t 

-0.10 

-0.20 

0.60 

0.00 

0.10  0.30 

0.05 

0.00 

Model  4 

-0.08 

-0.35 

0.22 

0.15 

Not  modeled 

— 

Model  5 

-0.06 

-0.46 

0.49 

0.29 

Not  modeled 

— 

Model  6 

-0.10 

-0.12 

0.63 

0.09 

Not  modeled 

— 

Model  7 

-0.08 

-0.46 

0.54 

0.26 

0.12  0.22 

0.03 

0.09 

Model  8 

-0.20 

-0.33 

0.02 

0.60 

-0.05  0.10 

0.64 

0.64 

Model  9 

-0.20 

-0.10 

0.62 

0.14 

0.03  0.28 

0.08 

0.08 

Model  9t 

-0.20 

-0.10 

0.65 

0.15 

0.04  0.29 

0.00 

0.08 

*depth  of  models  1  and  4  is  0.20m 


Table  23:  3D  Ellipsoid  inversion  results  (dimensions). 


60mm 

Error 

20mm 

Error 

A 

B 

C 

%Vol 

A  B 

C 

%Vol 

True 

0.051 

0.051 

0.051 

0.017  0.017 

0.017 

Model  1 

0.14 

0.08 

0.01 

2 

Not  modeled 

— 

Model  2 

0.13 

0.04 

0.03 

4 

Not  modeled 

— 

Model  3 

0.11 

0.04 

0.03 

1 

0.04  0.01 

0.01 

7 

Model  3t 

0.07 

0.05 

0.04 

1 

0.02  0.02 

0.02 

0 

True 

0.100 

0.030 

0.030 

0.035  0.010 

0.010 

Model  4 

0.10 

0.04 

0.03 

37 

Not  modeled 

— 

Model  5 

0.09 

0.06 

0.02 

351 

Not  modeled 

— 

Model  6 

0.21 

0.05 

0.03 

273 

Not  modeled 

— 

Model  7 

0.08 

0.04 

0.03 

12 

0.02  0.01 

0.01 

21 

Model  8 

0.01 

0.01 

0.01 

99 

0.14  0.05 

0.05 

8935 

Model  9 

0.12 

0.05 

0.05 

206 

0.05  0.02 

0.01 

191 

Model  9t 

0.13 

0.06 

0.04 

253 

0.05  0.01 

0.01 

61 
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The  best  results  were  obtained  using  the  3D  data  and  spherical  inversion  models.  The 
average  location  error  and  size  (radius)  error  were  approximately  half  those  of  the  2D 
data  with  spherical  models.  The  ellipsoid  inversion  required  a  greater  number  of  free 
parameters,  which  contributed  slightly  to  the  higher  errors  (Table  21  vs  22).  The  more 
significant  problem  is  believed  to  be  the  mis-match  of  ellipsoid  modeling  techniques. 
Model  9t  is  a  fully-populated  noise-free  version  of  Model  9.  Comparing  Models  9t  and 
3t  in  Tables  20-23  shows  that  even  with  perfect  data,  there  is  a  problem  with  matching 
ellipsoidal  models  that  is  not  present  for  spherical  models. 

It  was  also  found  that  the  raw  (uninterpolated)  data  points  produced  significantly  better 
results  than  the  gridded  results.  Comparable  tables  for  these  cases  are  presented  in 
Appendix  B-4.  The  interpolation  had  the  most  detrimental  effect  on  the  3D  data.  The 
original  300  data  points  were  interpolated  to  approximately  2,000  points,  many  of  which 
were  beyond  the  range  of  the  original  30  second  data  collection  window.  On  the  other 
hand,  the  2D  data  were  interpolated  from  205  to  400  points,  all  of  which  were  within  the 
outer  bounds  of  the  data  simulation.  There  was  also  a  difference  in  the  efficacy  of  the 
gridding  algorithm.  The  2D  Geosoft  minimum  curvature  routine  is  well  known  to 
produce  good  results  for  potential  fields.  The  performance  of  the  Rockware  algorithm  is 
relatively  untested.  A  comparison  between  the  two  was  conducted  and  the  Rockware 
gridding  produced  location  errors  between  3-5times  higher  than  the  Geosoft  grids  (see 
Appendix  B-4). 
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Conclusions 

In  all  cases,  the  3D  acquisition  and  processing  methodology  has  been  shown  to  be  better 
than  the  standard  2D  method  for  detection  and  location  of  closely  spaced  targets. 
Applying  1cm  of  positional  noise  and  O.lnT  of  magnetometer  noise  to  the  theoretical 
models,  the  best  techniques  averaged  location  errors  of  0.09m  and  depth  errors  of  0.06m. 
It  is  estimated  that  feedback  and  cooperative  inversion  techniques  can  further  reduce  this 
to  a  radial  error  of  0.05m. 

The  Gradient  String  approach  detected  all  of  the  targets  and  located  them  in  (x,  y)  more 
accurately  than  the  2D  analytic  signal  peaks.  Depth  estimates  using  the  analytic  signal 
ratio  and  Euler  deconvolution  were  five  times  more  accurate  using  the  3D  data.  The 
dipole  inversion  code  suffered  from  inadequate  signal  strength  and  segmentation 
problems.  This  caused  the  location  accuracies  to  fall  below  those  obtained  from  the 
simpler  methods.  The  3D  results,  however,  showed  significantly  fewer  problems  than  the 
2D  data.  The  solid-body  inversion  ran  more  smoothly  than  the  dipole  inversion,  but 
highlighted  a  potential  mis-match  in  theoretical  models  between  MAGMOD  and 
ModelVision.  The  best  results  were  produced  using  the  3D  data  and  spherical  models, 
where  target  locations  were  twice  as  accurate  as  the  comparable  2D  data  sets. 

It  was  also  shown  that  the  Rockware  3D  gridding  algorithm  was  not  as  well-suited  to 
potential  fields  as  the  Geosoft  approach.  Comparison  tests  showed  an  increase  in 
location  error  of  three  times.  The  3D  results  in  the  Detection,  AS  depth  and  Euler  depth 
analysis  sections  relied  on  the  Rockware  algorithm  and  were  as  good  or  better  than  the 
2D  data.  This  indicates  that  the  performance  enhancement  of  the  3D  method  is  still  being 
under  reported. 

This  study  did  not  examine  the  potential  benefits  of  the  3D  methodology  to  techniques 
for  discriminating  between  UXO  and  non-UXO.  This  was  due  to  the  limitations  of 
simulating  non-UXO  signatures.  Estimation  of  target  size  however,  a  key  element  in 
clutter  discrimination,  was  quite  accurately  detennined  for  the  spherical  target  sets. 


Transition 

This  project  has  found  considerable  advantage  to  a  3D  data  collection  and  processing 
methodology.  Transition  efforts  should  focus  on  (a)  added  benefits  than  can  be  obtained 
by  refinement  and  cooperative  use  of  the  techniques  explored  here  (b)  extending  the 
study  to  include  the  benefits  of  3D  methods  to  discrimination  of  UXO  from  non-UXO, 

(c)  introduction  of  commercial  software  and  hardware  vendors  to  improve  and  streamline 
the  techniques,  and  (d)  research  into  extension  of  these  techniques  to  electromagnetic 
data.  The  ultimate  objective  should  be  to  combine  these  techniques  into  a  hand-held 
sensor  platfonn  so  that  mapping,  detection  and  discrimination  can  be  conducted  in  a 
single  pass,  thereby  significantly  reducing  field  costs. 


37 


Recommendations 

There  are  several  additional  benefits  to  3D  processing  that  this  study  did  not  encompass. 
It  did  not  attempt  to  implement  any  sort  of  feedback  from  one  technique  to  another.  For 
example,  Euler  depth  estimates  can  be  used  to  refine  the  horizontal  location  and  then 
position  values  can  be  locked  to  reduce  the  number  of  free  parameters  in  the  inversion 
process.  This  would  free  up  variables  to  improve  the  resolution  of  target  shape  and  size. 
Similarly,  iterative  steps  through  the  inversion  to  fix  parameters  in  sequence  (invert  for 
location,  invert  for  size,  invert  for  shape,  invert  for  orientation)  is  another  process  that 
may  prove  useful.  The  use  of  irregularly  shaped  lattice  nodes  (e.g  10cm  x  10cm  x  5cm) 
may  offer  additional  benefits  by  reducing  the  buffer  zone  around  the  3D  lattice  so  that 
depth  calculations  can  be  made  a  lower  heights  with  better  signal-noise  ratios. 

This  study  should  be  extended  to  include  the  benefits  of  3D  methodology  to  current 
discrimination  techniques  (UXO  vs  non-UXO),  such  as  those  reported  by  University  of 
British  Columbia. 

Beyond  studies  into  the  benefits  of  3D  data  are  tests  of  the  practical  application  of  the 
method.  In  conjunction  with  the  collection  of  field  data  for  the  verification  phase  of  the 
project,  several  points  relating  to  software  need  to  be  addressed.  3D  gridding  routines 
need  to  be  improved  to  the  same  level  of  reliability  as  the  existing  2D  routines.  The 
interpolation  algorithms  in  the  existing  3D  software  were  shown  to  produce  locations 
errors  roughly  twice  that  of  the  2D  software.  The  extrapolation  of  data  beyond  the  3D 
“cloud”  in  order  to  fill  the  1x1x0. 5m  volume  highlighted  these  problems.  It  also 
indicates  that  additional  care  must  be  taken  to  collect  data  throughout  the  data  space  in 
order  to  maximize  the  utility  of  the  3D  method. 

The  discrepancies  between  ellipsoid  modeling  routines  needs  to  be  resolved  by 
examination  of  field  or  test-stand  data.  Cylindrical  ellipsoid  inversion  routines  must  also 
be  developed  to  overcome  some  of  the  software  shortcomings  discovered  in  this  study. 
Ultimately,  a  real-time  system  will  allow  technicians  to  digitally  map,  discriminate  and 
excavate  a  site  in  a  single  mobilization  without  repeated  referral  to  a  processing 
geophysicist.  This  would  considerably  reduce  the  time  and  effort  required  for  UXO 
investigations.  Development  of  any  new  software  routines  should  therefore  be  under  a 
single  processing  package  capable  of  hand-held  operation  for  the  field  staff,  and  desktop 
operation  for  QA  geophysicists.  This  will  also  facilitate  a  broader  distribution  and 
streamline  processing  during  the  initial  commercialization  period.  All  of  these  tasks  will 
require  the  involvement  of  relevant  commercial  software  and  hardware  vendors. 

The  final  task  is  to  extend  the  scope  of  this  study  to  include  electromagnetic  (EM) 
sensors.  This  will  be  sensor  specific,  probably  dealing  with  the  Geonics  EM-61HH  and 
the  Geophex  Gem-300.  Research  into  inversion  routines  and  the  effects  of  coil 
orientation  are  already  underway  through  SERDP/ESTCP  and  should  be  leveraged  with 
the  3D  methodology  proven  here.  This  will  allow  stand-alone  EM  techniques  to  be 
applied  in  noisy  magnetic  environments,  and  will  form  the  basis  for  future  joint  3D 
inversion  of  mag/EM  data. 
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Magnetic  modeling  and  inversion  routines  that  take  advantage  of  recorded  altitude  or  Z 
values  have  followed  the  development  of  improved  3 -dimensional  positioning 
technologies.  Usually,  this  technique  involves  recording  ground  clearance  (or  ellipsoid 
height  in  the  case  of  DGPS  systems)  during  the  course  of  collecting  routine  profile  data. 
Technologies  include  laser  or  radar  altimeters  for  airborne  techniques,  DGPS  heights,  or 
laser  and/or  ultra-sonic  positioning  for  ground  systems. 

Visualization  of  the  data,  however,  maintains  a  2D  approach  using  either  profiles  or 
grids.  Tools  such  as  3D  surface  maps  are  often  used  to  project  gridded  data  into  a 
simulated  3D  surface,  but  this  is  still  a  representation  of  2D  data. 

A  simple  approach  to  collecting  actual  3D  data  is  recording  multiple  sensors  at  different 
heights  during  a  single  profile,  such  as  vertical  magnetic  and  electromagnetic  gradients 
(Geometries  G-858  and  Geonics  EM-61).  The  utility  of  this  has  been  proven  in  many 
studies.  Many  inversion  routines  rely  on  gradients,  either  measured  or  calculated,  for 
depth  detennination  and  ordnance  discrimination. 

Alternatively  3D  data  may  be  acquired  with  a  single  sensor,  accurately  positioned,  and 
moved  around  in  a  volume  of  space  above  the  target  area  in  re-acquisition  mode.  This 
SERDP  funded  seed  project  investigates  this  approach  collection  of  data  in  a  3D  volume 
to  determine  the  resulting  improvement  in  inversion  capabilities. 
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Traditional  Approach 

The  traditional  approach  to  geophysical  mapping  for  UXO  has  been  to  collect  data  along  transects  at  a 
nominal  height  (1 D)  and  present  the  results  on  a  grid  (2D).  Single  profiles  show  the  data  with  great 
resolution  and  fidelity  since  they  are  not  subject  to  gridding  interpolation,  but  are  difficult  to  interpret.  The 
most  obvious  short-fall  being  the  lack  of  information  to  either  side  of  a  profile  anomaly.  By  presenting  the 
data  with  reference  to  neighboring  points  on  either  side,  the  2D  perspective  makes  it  much  simpler  to  pick 
out  discrete  targets  from  broader  anomalies.  This  simplification  is  an  important  part  of  the  interpretation 
process. 

Data  processing  techniques  based  on  2D  data  sets  also  have  significant  advantages  over  their  1 D 
counterparts.  There  are  many  advantages  to  analyzing  data  that  have  been  spatially  correlated  in  two 
dimensions.  Spatial  filters,  for  example,  can  remove  uncorrelated  noise  better.  FFT  processes,  such  as 
analytic  signal,  can  also  enhance  features  more  reliably  and  with  fewer  assumptions  when  applied  to 
grids  than  to  profiles. 

The  disadvantage  of  2D  presentations  is  that  it  assumes  that  all  of  the  data  were  collected  on  a  common 
plane,  although  this  is  not  usually  true.  Whenever  a  data  point  is  collected  at  different  height  than  the 
adjacent  data  points,  it  creates  an  artifact  in  the  grid.  Depending  on  the  severity  and  distribution  of  these 
artifacts,  any  interpretation  based  on  the  gridded  data  will  be  more  or  less  flawed. 

As  a  result,  many  additional  2D  data  processing  steps  must  be  incorporated  into  the  workflow  in  order  to 
correct  for  the  “common  plane”  assumption  and  produce  a  viable  interpretation.  In  addition  to  the 
gridding  process,  which  generally  desamples  data  down-line  and  interpolates  data  between  lines, 
filtering  or  deletion  of  data  points  must  be  applied  to  produce  a  satisfactory  grid.  In  some  cases,  these 
additional  steps  may  reduce  data  quality  more  than  they  enhance  it. 


State-of-the-Art  Approach 

The  impact  of  height  on  gridding  is  not  a  new  phenomenon.  The  effects  of  height  in  both  magnetic  and 
electromagnetic  data  are  well  documented.  Numerous  studies  have  been  done  to  either  limit,  quantify  or 
correct  these  effects.  Techniques  include  wheeled  carts  to  control  height,  altimeters  and  other 
positioning  systems  to  measure  height,  and  various  upward  or  downward  continuation  algorithms  to 
correct  readings  to  a  nominal  height. 

Advances  in  positioning  systems  and  computing  capabilities  now  allow  collection  and  storage  of 
geophysical,  position  and  height  data  in  a  single  database.  Software  packages  for  data  inversion  are 
taking  advantage  of  increased  power  and  are  inverting  data  with  height  information.  The  MTADS 
inversion  code  and  the  commercially  available  ModelVision  from  Encom  are  two  examples. 

Although  the  technology  has  incorporated  several  aspects  of  3D  data,  all  of  these  are  attempts  to  solve 
the  problem  of  continually  forcing  3D  data  into  a  2D  format.  To  date,  the  best  inversion  solutions  ignore 
the  2D  format  altogether  and  work  entirely  with  profile  data.  Although  this  work  flow  avoids  many  of  the 
problems,  it  does  not  take  advantage  of  any  of  the  positive  elements  of  2D  processing  such  as  noise 
reduction  or  FFT  algorithms.  It  makes  use  of  only  one  dimension  of  spatial  correlation  (down-line).  Only 
filters  and  processes  that  can  be  applied  to  profiles  are  used. 

Notes  on  Dimensions 

For  descriptive  clarity,  a  distinction  must  be  made  between  the  number  of  dimensions  of  position  data 
collected  during  a  survey  and  the  number  of  dimensions  of  the  presentation  format.  The  traditional 
survey  measures  two  dimensions  of  position  data  (X,Y)  and  presents  the  geophysical  results  as  1 D 
profiles  (a  projection  against  a  single  axis  of  time  or  horizontal  position,  with  the  data  as  the  second  axis), 
2D  grids/contours  (a  projection  against  two  horizontal  position  axes  using  color  to  represent  the  data 
axis),  or  a  3D  surfaces  /curves  / nets  (a  projection  against  two  horizontal  position  axes,  using  the  illusion  of 
a  third  vertical  axis  to  represent  the  data).  Data  may  also  be  presented  as  profiles  along  survey  lines 
plotted  on  2D  maps,  which  are  similarto  the  2D  grids. 

Even  with  the  recording  of  3D  positional  data  (X,Y  and  height),  none  of  the  aforementioned  presentation 
techniques  display  or  accommodate  the  additional  height  information.  This  is  a  limitation  of  all  two- 
dimensional  media  (paper  and  screens).  It  is  extremely  difficult  to  represent  3D  positions  on  a  2D  screen 
or  map.  This  does  not  include  a  dimension  for  the  geophysical  data  itself.  Even  3D  nets  are  technically 
2D  presentations  using  illusions  of  perspective  and  shading  to  appear  3D.  Occasionally,  a  3D 
presentation  is  rotated  to  enhance  the  illusion,  thus  adding  time  as  an  added  dimension  to  a  2D  image. 
The  3D  holographic  presentations  of  science  fiction  would  still  require  rotation  of  the  image  or  movement 
of  the  observer  in  order  to  fully  interpret  it,  just  as  any  solid  object  would. 

Multi-dimensional  data  sets  therefore  require  that  we  either  reduce  the  information  to  it’s  most  relevant 
form  for  presentation,  or  else  to  think  and  to  present  data  in  terms  of  matrices  and  functions  with  n- 
dependent  variables.  This  applies  equally  to  3D  positioned  data  from  a  single  sensor,  multiple  responses 
from  a  single  sensor  (ie  EM  time  gates),  or  integrated  multi-sensor  data  sets.  Although  this  makes  them 
ideal  candidates  for  computer  analysis,  which  is  based  almost  entirely  on  linear  matrices,  the  human 
factor  is  still  required  for  control  and  judgement. 
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Background 

This  projectfollows  from  the  following  points: 

•  traditional  survey  techniques  lose  much  of  the  benefit  of  2D  data  processing  technic 
require  additional  processing  to  remove  or  diminish  the  effects  of  having  forced  o 
plane 

•  state-of-the-art  processing  techniques  that  incorporate  3D  positional  information 
benefit  of  2D  data  processing  by  ignoring  the  potential  of  spatially  correlating  the 
straight  to  numerical  inversion  after  1 D  processing 

•  data  collected  with  3D  positional  information  offers  the  potential  of  spatially  co 
vertically  as  well  as  horizontally 

This  last  point  presents  particularly  interesting  possibilities.  It  deals  with  the  problerr 
changes  in  height  by  handling  3D  data  in  three-dimensions,  as  well  as  providing  opport 
the  state-of-the-art  inversion  techniques  by  incorporating  the  benefits  of  spatial  c 
inversion. 


Approach 

For  developmental  purposes,  the  approach  to  creating  and  analyzing  these  new  tech 
the  data  reacquisition  process.  After  an  anomaly  of  interest  is  selected,  reacquisition 
with  a  similar  instrument  using  a  3D  positioning  system.  The  effect  is  to  collect  a  “clc 
over  a  suspect  anomaly. 

The  “cloud”  of  data  are  then  interpolated  into  a  3D  lattice  (grid)  for  processing.  T 
opportunity  for  examination  of  3D  analogs  of  2D  processing  techniques.  Several  diffi 
techniques  are  available.  Two  of  the  most  accessible  software  packages  forthis  purpc 
RockWorks.  Matlab,  however,  only  offers  linear  and  nearest-point  lattice  routines, 
eight  methods,  six  of  which  are  suitable  forthis  application. 

3D  Filters  and  Products 

3D  data  processing  techniques  can  now  be  applied  to  the  data  prior  to  analysis.  These 
analogs  to  standard  2D  techniques,  such  as  windowing  and  filtering,  as  well  as  me 
techniques  to  derived  products  such  as  vertical  gradient  and  analytic  signal.  The  ac 
calculations  of  gradient  products,  of  course,  is  the  fact  that  far  more  accurate  vertical  g 
hand  than  ever  before.  Higher  order  gradients  are  also  possible  as  directly  measured  f 
than  mathematical  extrapolations  of  total  field  data.  New  filters  and  new  products  are 
lattice  data.  Moving  average  filters  with  variable  weighting  based  on  absolute  or  rela 
used  to  take  advantage  of  the  variable  signal-to-noise  ratios  produced  by  ground  cluttei 

3D  Gradient  Strings 

Although  inversion  routines  do  not  rely  on  visualization  of  the  data,  reduction  of  the 
intrinsic  qualities  is  essential  to  presentation  for  human  interpretation  and  control  of  filt< 
inversion  parameters.  Several  standard  presentation  techniques  are  reviewed,  a 
concept  of  Gradient  Strings  is  introduced. 

Standard  presentation  techniques  of  lattice  data  are  layers,  slices,  iso-surfaces  an 
Gradient  strings,  lines  through  the  lattice  tracing  maximum  rate  of  change  in  the  field 
introduced  here.  It  is  applicable  to  both  magnetic  and  electromagnetic  data,  j 
interpretation  of  the  physics  differs,  the  utility  for  processing  and  target  identification  is  tl 


The  visualization  and  manipulation  of  gradient  strings  will  allow  the  operator  to 
parameters  and  to  deselect  targets  from  noise  or  to  pin-point  multiple  targets  within  clc 
underlying  assumption  is  that  maximum  gradient  lines  for  legitimate  sources  will  be  o 
increase  in  magnitude  with  proximity  to  the  ground.  Noise,  which  might  otherwise  be  [ 
targets,  will  lack  this  continuity  through  3D  space.  Similarly,  multiple  targets  that  may ; 
target  at  survey  heights,  can  be  resolved  as  distinct  gradient  strings.  Processing 
gradient  string  information  can  include  spatial  windows,  cut-off  thresholds,  and  filter  str 
of  peaks  can  be  controlled  by  length  and  continuity  of  gradient  strings. 
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Analysis  of  the  benefits  of  3D  processing  will  be  based  on  a  suite  of  magnetic  dipole  models.  Although 
this  technique  is  equally  applicable  to  electromagnetic  models,  magnetic  modeling  offers  fewer  variables 
and  more  robust  inversion  at  this  stage.  Models  include  single  dipoles  at  multiple  depths,  and  pairs  of 
dipoles  of  varying  strengths,  depths  and  offsets.  A  duplicate  set  of  models  for  filter  testing,  with  random 
noise  added  to  the  positions  and  geophysical  readings  has  also  been  created.  Having  developed  the 
theoretical  models,  1 D  profiles  and  2D  grids  are  extracted  and  inverted  using  different  techniques  to 
establish  a  baseline  performance  for  each  of  the  inversion  methods.  These  methods  will  include 
amplitude  thresholds,  Euler  deconvolution,  the  MTADS  and  ModelVision  inversion  codes. 

3D  Gridding 

Various  3D  gridding  algorithms  are  being  tested  to  determine  which  produced  the  most  acceptable 
results.  Results  at  lattice  nodes  will  be  comapred  to  the  theoretical  response  to  determine  the  quality  of  a 
given  routine.  Input  data  will  begin  with  a  fully  sampled  space.  Data  density  will  then  be  reduced  in 
stages.  Noise  can  then  be  added  to  provide  a  targetforfiltering  routines. 

3D  Processing 

3D  analogs  to  2D  processing  procedures  are  then  developed.  These  have  included  Hanning  filters  and 
Analytic  Signal  calculations.  More  advanced  algorithms  for  weighted  filters  and  gradient  string 
calculations  have  also  been  developed.  Modifications  to  the  various  presentation  formats  in  different 
software  suites  are  being  explored  in  order  to  streamline  analysis  and  itemize  features  for  future 
development. 


Inversion 

The  various  inversion  methods  are  then  tested  against  the  synthetic  3D  data  (with  and  without  noise)  and 
compared  to  the  baseline  results  to  quantify  any  improvement.  This  will  require  some  modification  to  the 
standard  procedures  for  some  of  the  techniques.  The  amplitude  threshold  is  the  simplest  method.  In 
three  dimensions,  it  is  an  extension  of  the  gradient  string  principle.  Rather  than  picking  a  peak  from  a 
grid,  the  peak  is  traced  through  3D  space  and  can  include  a  cut-off  for  length  as  well  as  amplitude.  The 
Euler  routine  will  include  the  measured  gradient  rather  than  a  calculated  one  in  its  algorithm.  The  MTADS 
and  ModelVision  codes  operate  on  profile  lines.  These  will  be  extracted  from  the  3D  lattice  and  included 
in  the  inversion  database. 
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Gradient  String  Demo 

This  is  a  simple  demonstration  of  the  gradient  string 
concept.  Analytic  Signal  data  are  presented  as  a 
transparency,  with  peaks  shown  as  solids  at  the 
appropriate  node.  The  line  of  peaks  at  the  center  of  the 
lattice  are  continuous  and  increase  in  intensity  with 
decreasing  height.  The  single  point  peaks  are  the 
result  of  noise  in  the  lattice  data.  Gridded  data 
(horizontal  slices)  would  intersect  one  or  more  local 
lattice  peaks  and  might  be  misinterpreted.  Using  the 
3D  spatial  relationship  between  points  allows  for  less 
ambiguous  interpretation.  These  peaks  can  be  used  to 
set  filter  thresholds  in  processing,  amplitude  thresholds 
in  peak  selection,  or  continuity  thresholds  in  gradient 
string  selection. 
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Data  Presentation 

At  right  is  a  screen  shot  from  ModelVision 
demonstrating  3D  inversion  of  synthetic  magnetic 
data.  2D  gridded  data  in  the  top  left  shows 
magnetic  response  from  lowest  layer  of  lattice. 
Superimposed  are  two  spheres  presented  as  the 
starting  model  of  the  inversion  process. 
Parameters  of  the  spheres  are  shown  in  the 
bottom  right. 

The  software  is  semi-manual,  in  that  the  user 
must  provide  a  reasonable  starting  model  before 
proceeding  to  inversion.  Inversion  allows  for  free 
parameters  such  as  position,  depth,  radius, 
susceptibility  and  remanence  of  spheres. 
Additional  shapes,  such  as  spheroids  and  various 
regular  and  irregular  polygons  are  supported  with 
variations  in  free  parameters  as  required. 

Below  the  grid  are  two  NS  profiles,  also  from  the 
lowest  level  of  the  lattice.  Each  one  passes  over 
one  of  the  spheres.  The  black  profile  is  the 
original  data,  the  red  profile  is  the  starting  model 
for  the  inversion  process. 

The  3D  representation  in  the  upper  right  shows 
the  two  spheres  and  the  original  total  field  data 
corresponding  to  the  grid  and  profiles  on  the  left. 
Line  paths  representing  three  layers  of  the 
complete  lattice  are  also  shown. 
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Abstract 

Magnetic  modeling  and  inversion  routines  that  take  advantage  of  recorded  altitude  or  Z 
values  have  followed  up  on  the  development  of  improved  3-dimensional  positioning 
technologies.  Usually,  this  technique  involves  recording  ground  clearance  or  ellipsoid  height 
during  the  course  of  collecting  routine  profile  data.  Technologies  include  laser  or  radar 
altimeters,  DGPS  heights,  or  laser  and/or  ultra-sonic  positioning  systems.  Visualization  of  the 
data,  however,  maintains  a  2D  approach  using  either  profiles  or  grids.  This  is  an  obvious 
limitation  of  2D  media  such  as  paper  and  computer  screens.  Tools  such  as  3D  surface  maps  are 
often  used  to  project  gridded  data  into  simulated  3D,  but  this  is  not  the  same  thing  as  presenting 
3D  data.  Processing  of  these  data  often  takes  a  step  backwards  and  is  limited  to  ID  tools.  The 
noise  reduction  and  interpretive  capabilities  derived  from  spatially  correlating  data  across  lines  is 
lost. 

This  paper  presents  a  justification  for  3D  collection,  processing  and  interpretation  of  data 
for  small,  discrete  objects  such  as  drums,  pipes  and  unexploded  ordnance  that  are  common  to 
many  environmental  surveys.  In  addition  to  traditional  visualization  techniques,  the  concept  of 
“gradient  strings”  and  their  utility  in  data  analysis  is  introduced.  While  this  paper  focuses  on 
magnetic  data,  the  basic  principles  are  equally  applicable  to  electromagnetic  data. 

Introduction 

Many  geophysical  data  sets,  collection  procedures  and  presentation  techniques  are 
described  as  3-dimensional.  For  the  purposes  of  clarity,  a  distinction  must  be  made  between  the 
number  of  dimensions  of  the  position  data  collected  during  the  survey  and  the  number  of 
dimensions  of  the  presentation  format.  The  traditional  survey  measures  two  dimensions  of 
position  data  (X,Y)  and  presents  the  geophysical  results  as: 

•  ID  profiles  (a  projection  against  a  single  axis  of  time  or  horizontal  position,  with  the 
data  as  the  second  axis), 

•  2D  grids/contours  (a  projection  against  two  horizontal  position  axes  using  color  to 
represent  the  data  axis), 

•  or  3D  surfaces  /curves  /nets  (a  projection  against  two  horizontal  position  axes,  using 
the  illusion  of  a  third  vertical  axis  to  represent  the  data). 

Data  may  also  be  presented  as  profiles  along  survey  lines  plotted  on  2D  maps,  which  are 
similar  to  the  2D  grids,  although  this  format  has  apparently  fallen  out  of  fashion  since  the  advent 
of  simple  gridding  and  contouring  packages. 

Even  with  the  recording  of  3D  positional  data  (X,Y  and  height),  none  of  the 
aforementioned  presentation  techniques  display  or  accommodate  the  additional  height 
information.  Generally,  the  height  is  assumed  to  be  constant  and  the  data  are  forced  into  a  2D 
grid.  This  is  a  limitation  of  all  two-dimensional  media  (paper  and  screens).  Even  3D  surfaces 
are  technically  2D  presentations  using  illusions  of  perspective  and  shading  to  appear  3D.  It  is 


extremely  difficult  to  represent  3D  points  on  a  2D  screen  or  map,  and  this  does  not  yet  include  a 
dimension  for  the  geophysical  data  itself.  Occasionally,  the  perspective  of  a  3D  presentation  is 
rotated,  thus  adding  time  to  enhance  the  illusion  of  dimension.  Exchanging  our  2D  media  for  the 
3D  holographic  presentations  of  science  fiction  would  still  require  some  form  of  transparency 
plus  rotation  of  the  image  or  movement  of  the  observer  in  order  to  fully  interpret  it,  just  as  any 
solid  object  would.  A  3-dimensional  representation  therefore  presents  several  fundamental 
obstacles  to  the  observer. 

Multi-dimensional  data  sets  therefore  require  that  we  either  reduce  the  infonnation  to  its 
most  relevant  form  for  presentation,  or  else  to  think  and  to  present  data  in  terms  of  matrices  and 
functions  with  n-dependent  variables.  This  applies  equally  to  3D  positioned  data  from  a  single 
sensor,  multiple  responses  from  a  single  sensor  (ie  EM  time  gates),  or  integrated  multi-sensor 
data  sets.  Although  this  makes  them  ideal  candidates  for  computer  analysis,  which  is  based 
almost  entirely  on  linear  matrices,  the  human  factor  is  still  required  for  control  and  judgment. 

Traditional  Data  Collection 

A  single  reading  at  one  point  in  space  rarely  has  any  geophysical  significance.  Until  the 
reading  is  presented  with  respect  to  neighboring  points  it  has  little  relevance.  It  is  this  spatial 
correlation  that  gives  a  data  profile  meaning.  Similarly,  consecutive  profiles  allow  the  data  to  be 
spatially  correlated  in  two  dimensions  -  down  line  and  across  line.  The  traditional  approach  to 
geophysical  mapping  therefore  is  to  collect  data  along  transects  (ID)  at  a  nominal  height  and 
present  the  results  as  a  grid  (2D).  Each  added  dimension  brings  a  new  level  of  relevance  to  the 
interpretation,  but  there  are  drawbacks. 

Single  profiles  show  the  data  with  great  resolution  and  fidelity  since  they  are  not  subject 
to  gridding  interpolation,  but  are  difficult  to  interpret.  The  most  obvious  short-fall  is  the  lack  of 
information  to  either  side  of  a  profile  anomaly.  By  presenting  the  data  with  reference  to 
neighboring  points  on  either  side,  the  2D  perspective  makes  it  much  simpler  to  pick  out  discrete 
targets  from  broader  anomalies.  This  simplification  is  an  important  part  of  the  interpretation 
process. 

Data  processing  techniques  based  on  2D  data  sets  also  have  significant  advantages  over 
their  ID  counterparts.  There  are  many  advantages  to  analyzing  data  that  have  been  spatially 
correlated  in  two  dimensions.  Spatial  filters,  for  example,  can  remove  uncorrelated  noise  better. 
FFT  processes,  such  as  analytic  signal,  can  also  enhance  features  more  reliably  and  with  fewer 
assumptions  when  applied  to  grids  than  to  profiles. 

The  disadvantage  of  2D  presentations  is  that  it  assumes  that  all  of  the  data  were  collected 
on  a  common  plane,  although  this  is  not  usually  true.  Whenever  a  data  point  is  collected  at 
different  height  than  the  adjacent  data  points,  it  creates  an  artifact  in  the  grid.  Depending  on  the 
severity  and  distribution  of  these  artifacts,  any  interpretation  based  on  the  gridded  data  will  be 
more  or  less  flawed. 

As  a  result,  many  additional  2D  data  processing  steps  must  be  incorporated  into  the  work 
flow  in  order  to  correct  for  the  “common  plane”  assumption  and  produce  a  viable  interpretation. 
In  addition  to  the  gridding  process,  which  generally  desamples  data  down-line  and  interpolates 
data  between  lines,  filtering  or  deletion  of  data  points  must  be  applied  to  produce  a  satisfactory 
grid.  In  some  cases,  these  additional  steps  may  reduce  data  quality  more  than  they  enhance  it. 


State-of-the-Art  Data  Collection 


The  impact  of  height  on  gridding  is  not  a  new  phenomenon.  The  effects  of  height  in  both 
magnetic  and  electromagnetic  data  are  well  documented.  Numerous  studies  have  been  done  to 
limit,  quantify  or  correct  these  effects.  Techniques  include  wheeled  carts  to  control  height, 
altimeters  and  other  positioning  systems  to  measure  height,  and  various  upward  or  downward 
continuation  algorithms  to  correct  readings  to  a  nominal  height. 

Advances  in  positioning  systems  and  computing  capabilities  now  allow  collection  and 
storage  of  geophysical,  position  and  height  data  in  a  single  database.  Software  packages  for  data 
inversion  are  taking  advantage  of  increased  power  and  are  inverting  data  with  height 
information. 

Although  the  technology  has  incorporated  several  aspects  of  3D  data,  all  of  these  are 
attempts  to  solve  the  problem  of  continually  forcing  3D  data  into  a  2D  format.  To  date,  the  best 
inversion  solutions  ignore  the  2D  format  altogether  and  work  entirely  with  profile  data. 
Although  this  work  flow  avoids  many  of  the  problems,  it  does  not  take  advantage  of  any  of  the 
positive  elements  of  2D  processing  such  as  noise  reduction  or  FFT  algorithms.  It  makes  use  of 
only  one  dimension  of  spatial  correlation  (down-line).  Only  filters  and  processes  that  can  be 
applied  to  profiles  are  used. 


3D  Data  Collection 

The  only  way  to  take  advantage  of  spatial  correlation  of  data  collected  with  3D 
positioning,  is  to  process  it  in  three  dimensions.  Not  only  does  this  effectively  deal  with  the 
problems  associated  with  changes  in  height,  but  it  offers  the  potential  of  spatially  correlating  the 
data  vertically  as  well  as  horizontally.  This  also  provides  an  opportunity  to  enhance  the  state-of- 
the-art  inversion  techniques  by  incorporating  the  benefits  of  spatial  correlation  prior  to  inversion. 

For  logistical  and  testing  purposes,  this  paper  assumes  that  the  data  are  collected  as  part 
of  a  cued  investigation  (to  follow-up  a  previously  detected  anomaly)  using  magnetometer 
sensors,  although  the  techniques  are  equally  applicable  to  traditional  survey  collection  and  to 
electromagnetic  data  after  additional  experimentation  with  3D  gridding  and  inversion  techniques. 
After  an  initial  geophysical  survey,  specific  anomalies  are  identified  for  additional  follow-up. 
Prior  to  excavation,  the  operator  collects  a  “cloud”  of  data  points  at  irregularly  spaced  intervals 
around  the  pre-detennined  location.  The  processing  and  inversion  are  then  applied  to  this  data 
set  (Figure  1). 

The  “cloud”  of  data  is  then  interpolated  into  a  3D  lattice  (grid)  for  processing.  This 
offers  the  first  opportunity  for  examination  of  3D  analogs  of  2D  processing  techniques.  Several 
different  interpolation  techniques  are  available.  Two  of  the  most  accessible  software  packages 
for  this  purpose  are  Matlab  and  RockWorks.  Matlab,  however,  only  offers  linear  and  nearest- 
point  lattice  routines.  Rockworks  offers  eight  methods,  six  of  which  are  suitable  for  this 
application.  The  most  reliable  of  these  was  their  “inverse  distance/anisotropic”  technique.  This 
method  interpolates  nodes  based  on  an  average  of  the  closest  points  within  each  octant  weighted 
by  their  proximity  to  the  interpolated  node  (Figure  2). 
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Figure  1:  Total  field  magnetic  data  (in  nT)  collected  as  a  “cloud”  of  points  in  three  dimensions 
above  two  synthetic  targets.  Data  are  presented  as  a  10cm  cube  at  each  reading  location.  Data 
space  is  300cm  x  300cm  x  60cm.  Sample  density  is  l/6th  of  the  total  volume  (comparable  to  a 
fully  populated  2D  grid  over  the  same  area). 


The  targets  for  the  sample  data  set  presented  here  are  two  spherical  dipoles  with  different 
sizes  and  depths  in  close  proximity  to  one  another.  The  large,  deep  target  appears  in  the  SW 
quadrant  of  the  figures,  while  the  smaller,  shallower  target  appears  in  the  NE  quadrant. 
Additional  target  details  are  provided  in  Table  1.  The  collected  data  set  ranges  from  -100cm  to 
+200cm  in  X  and  Y,  and  from  10cm  to  60cm  above  ground  level.  It  will  be  seen  that  the  sensor 
height  is  a  critical  element  of  the  interpretation.  The  response  of  the  larger  target  is  visible  at  all 
heights,  but  the  rapid  increase  in  gradient  of  the  shallower  target  makes  it  much  more  definitive 
at  lower  heights.  Since  cued  investigations  typically  use  analog  instruments  held  very  close  to 
the  ground,  it  is  a  debatable  and  subjective  call  as  to  which  target  would  actually  be  excavated  by 
a  field  crew  -  the  larger  target,  the  stronger  response,  or  both.  Unless  both  are  excavated  and 
documented,  the  measured  data  will  be  difficult  to  use  as  a  calibration  tool  to  judge  system 
performance,  future  results  and/or  QA  assessment. 
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Figure  2:  Total  field  magnetic  data  (nT)  interpolated  from  “cloud”  of  points  in  figure  1  using  an 
inverse  distance  algorithm.  Data  space  is  300cm  x  300cm  x  60cm. 


This  3D  data  set  has  been  sampled  at  random  locations  throughout  this  data  space  at  l/6th 
the  total  volume  so  as  to  be  more  directly  comparable  to  a  fully  populated  2D  grid  over  the  same 
area.  At  this  density  and  distribution,  the  3D  gridding  routine  demonstrates  a  0.32nT  rms  error 
between  the  interpolated  results  and  the  true  modeled  values.  Gaussian  noise  of  2.5cm  was 
added  to  the  XY  positions,  5.0cm  to  the  Z  positions,  and  O.lnT  was  added  to  the  magnetic  data 
prior  to  3D  interpolation  and  processing. 

3D  data  processing  techniques  can  now  be  applied  to  the  data  prior  to  analysis.  These 
include  simple  3D  analogs  to  standard  2D  techniques,  such  as  windowing  and  filtering,  as  well 
as  modified  processing  techniques  to  derived  products  such  as  vertical  gradient  and  analytic 
signal.  For  example,  the  coefficients  for  a  3-point  symmetric  3D  Hanning  filter  are  provided  in 
Table  2.  The  advantage  to  lattice  calculations  of  gradient  products,  of  course,  is  the  fact  that  far 
more  accurate  vertical  gradient  data  are  at  hand  than  ever  before.  Higher  order  gradients  are  also 
possible  as  directly  measured  parameters,  rather  than  mathematical  extrapolations  of  total  field 
data.  New  filters  and  new  products  are  also  possible  with  lattice  data.  Moving  average  filters 
with  variable  weighting  based  on  absolute  or  relative  height  can  be  used  to  take  advantage  of  the 
variable  signal-to-noise  ratios  produced  by  ground  clutter.  The  result  is  a  data  set  that  has  been 
enhanced  with  spatial  correlation  in  all  three  directions  without  requiring  excessive  filters  to 
correct  for  2D  assumptions. 


Table  1:  Synthetic  target  parameters  used  for  sample  data  set.  Targets  were  modeled  as 
spherical  dipoles  with  a  magnetization  of  120  A/m  in  a  field  of  50,000nT  at  60°  inclination  and 
20°  West  declination. 


ID 

X  (cm) 

Y  (cm) 

Z  (cm) 

radius  (cm) 

1 

0.0 

0.0 

100.0 

5.0 

2 

100.0 

100.0 

20.0 

1.5 

Table  2:  Coefficients  of  3-point  3D  Hanning  filter 
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Figure  3:  Total  field  magnetic  data  (nT)  after  application  of  a  3 -point,  3 -dimensional  Hanning 
filter.  Data  space  is  300cm  x  300cm  x  60cm.  (top-left)  solid  volume  shows  faces  only,  (top- 
right)  horizontal  slice  at  40cm  above  the  ground,  (bottom-left)  two  pairs  of  iso-surfaces  showing 
volumes  below  50,000nT,  and  above  50,002nT.  (bottom-right)  the  same  iso-surfaces  as  before, 
with  a  75%  transparency. 


Visualization  techniques  must  now  be  applied  to  view  the  interior  of  the  3D  data  set. 
Several  different  techniques  are  demonstrated  in  Figure  3.  These  include  solids,  slices,  iso¬ 
surfaces  and  transparencies.  Slices  are  simply  planes  cut  through  the  lattice  at  varying  angles. 
Multiple  vertical  slices  are  often  called  fences.  Horizontal  slices  are  equivalent  to  2D  grids, 
except  that  they  can  have  been  interpolated  onto  a  plane,  rather  than  assumed  to  have  been 
collected  that  way.  Iso-surfaces  are  the  3D  equivalent  of  2D  contour  lines.  They  appear  as 
three-dimensional  blobs  encompassing  all  values  greater  than  the  specified  threshold.  Unlike 
contour  maps  only  a  single  surface  value  can  be  displayed  at  a  time,  although  multiple  iso¬ 
surfaces,  such  as  the  positive/negative  lobes  demonstrated  in  Figure  3,  are  possible.  Successive 
increments,  fluctuations  or  depressions  within  the  iso-surfaces  cannot  easily  be  represented. 
Transparencies  attempt  to  do  this  by  taking  the  solid  colors  and  reduce  their  intensity,  allowing 
visual  space  for  underlying  colors  and  shapes  to  be  presented.  All  of  these  techniques  require 
interactive  manipulation  in  order  to  be  useful  to  the  interpreter.  This  requires  additional  time  and 
effort  that  is  ill-suited  to  a  production  setting  where  hundreds  of  targets  must  be  analyzed. 

3D  Gradient  Strings 

Although  inversion  routines  do  not  rely  on  visualization  of  the  data,  reduction  of  the  data 
to  their  most  intrinsic  qualities  is  essential  to  presentation  for  human  interpretation  and  control  of 
filter,  processing  and  inversion  parameters.  Gradient  strings  are  a  new  concept  introduced  here 
to  simplify  the  visualization  of  3D  data  sets  typical  in  magnetic  and  electromagnetic  detection  of 
UXO  and  other  discrete  environmental  targets.  Although  the  interpretation  of  the  physics  differs 
between  magnetics  and  electromagnetics,  the  utility  for  processing  and  target  identification  is  the 
same.  These  are  essentially  lines  through  the  lattice  tracing  maximum  rate  of  change  in  the  field. 
For  magnetic  data,  one  can  imagine  layers  of  analytic  signal  grids  with  a  line  connecting  the 
peaks  from  one  layer  to  the  next. 

The  aim  of  this  technique  is  to  reduce  the  data  set  to  its  most  important  qualities  and 
allow  simple  visualization  with  as  few  control  parameters  as  possible.  The  underlying 
assumption  is  that  maximum  gradient  lines  for  legitimate  sources  will  be  continuous  and  will 
increase  in  magnitude  with  proximity  to  the  ground.  Noise,  which  might  otherwise  be  picked  as 
a  possible  target,  will  lack  this  continuity  through  3D  space.  Similarly,  multiple  targets  that  may 
appear  as  a  single  target  on  the  original  survey  data  can  be  resolved  as  distinct  gradient  strings. 
Processing  options  based  on  gradient  string  information  can  include  spatial  windows,  cut-off 
thresholds,  and  filter  strengths.  Selection  of  peaks  can  be  controlled  by  length  and  continuity  of 
gradient  strings.  If  the  line  branches  with  depth,  the  operator  is  faced  with  the  possibility  of 
multiple  targets  or  a  single  large  target  in  the  very  near  surface. 

For  a  string  to  be  valid,  therefore,  it  must  (a)  exceed  an  intensity  threshold  (b)  be 
continuous  over  several  layers  and  (c)  increase  in  intensity  with  depth.  The  result  is  a  line 
through  the  3D  data  space  that  points  directly  to  the  target.  To  further  aid  in  interpretation,  the 
line  is  color-coded  to  gradient  intensity. 

These  gradient  strings  and  the  numerical  inversion  provide  mutual  feedback  into  the 
decision  making  process.  The  visualization  and  manipulation  of  gradient  strings  will  allow  the 
operator  to  refine  processing  parameters  and  to  deselect  targets  from  noise  or  to  pin-point 
multiple  targets  within  close  proximity  to  one  another.  In  return,  the  inversion  will  provide 
quantitative  results  in  the  form  of  fit  coefficients. 
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Figure  4:  Analytic  signal  (nT/cm)  calculated  from  the  3D  total  field  data.  Slices  represent  data 
at  10cm,  30cm  and  60cm  above  the  ground  level.  Data  space  is  300cm  x  300cm  with  an 
expanded  vertical  scale.  Note  the  effect  of  second  order  gradient  effects.  The  smaller 
(shallower)  target  in  the  NE  quadrant  shows  a  much  stronger  gradient  response  at  lower  heights, 
but  has  almost  disappeared  at  greater  heights,  although  the  larger  (deeper)  target  is  visible  at  all 
heights. 
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Figure  5:  Gradient  strings  of  analytic  signal  (nT/cm).  (top-left)  all  strings  regardless  of  length, 
amplitude  or  continuity,  (top-right)  strings  after  application  of  0.3nT/cm  threshold,  (bottom- 
left)  strings  with  three  or  more  points,  regardless  of  amplitude,  (bottom-right)  strings  with  three 
or  more  points  and  amplitudes  greater  than  0.3nT/cm. 


3D  Inversion 


To  compare  inversion  results  from  2D  vs  3D  data  collection,  three  different  data 
collection  scenarios  were  modeled  on  the  sample  target  set  used  above.  The  first  was  a  set  of 
parallel  survey  lines  at  0.5m  spacing  and  0.4m  height,  the  second  was  a  2D  grid  at  0.1m  spacing 
and  0.4m  height,  and  the  third  was  a  3D  lattice  at  randomly  sampled  locations.  The  number  of 
data  points  for  the  2D  grid  and  the  3D  lattice  were  the  same.  Each  scenario  was  inverted  using 
the  Encom  ModelVision  software.  The  same  poorly  fit  starting  model  was  used  for  all 
inversions.  Results  of  the  inversion  are  presented  in  Table  3.  Figure  6  shows  a  screen  capture  of 
the  inversion  software  for  the  noise-free  case. 

Adding  a  small  amount  of  noise  (Gaussian  noise  2.5cm  to  XY,  5.0cm  to  Z,  O.lnT  to 
magnetics),  produced  correspondingly  poorer  results.  If  the  0.5m  survey  lines  crossed  directly 
over  the  two  targets,  the  inversion  fit  was  2.8%  rms.  If  the  survey  lines  are  offset  by  10cm  from 
the  target  center,  the  error  increases  to  3.6%  rms.  Given  that  most  geophysical  surveys  use  lm 
line  spacing  (rather  than  the  0.5m  used  here)  and  experience  2-5x  the  positional  error,  this  is 
probably  better  than  most  traditional  survey  results  can  produce.  The  fully  populated  2D  grid 
inverted  with  a  fit  of  2.9%  error,  while  the  3D  lattice  inverted  with  a  fit  of  2.2%  error. 

In  each  case,  the  3D  lattice  produced  slightly  better  numerical  inversion  results  than  the 
fully  populated  2D  grid  using  the  same  number  of  data  points.  This  reflects  the  predictive  nature 
of  potential  fields.  An  electromagnetic  system  would  presumably  show  an  even  greater 
difference  in  performance.  It  must  be  recognized,  however,  that  this  is  not  a  typical  2D  magnetic 
grid.  It  is  a  fully  populated  grid  with  a  raw  data  spacing  of  10cm  and  tight  vertical  controls. 
Since  this  is  best  accomplished  in  practice  using  a  3D  positioning  system  and  the  data  collection 
procedure  outlined  above,  the  true  measure  of  the  3D  inversion  results  for  field  applications  is 
the  comparison  to  the  10cm  offset  lines. 


Table  3:  Inversion  parameters  and  resul 


ts.  The  same  starting  model  was  used  for  all  inversions. 
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True 

0 

0 

100 

5.0 

0.24 

100 

100 

20 

1.5 

0.24 

0.02 

Start 

10 

-10 

20 

2.0 

0.30 

90 

110 

20 

2.0 

0.30 

18.4 

Noise-free 

Lines* 

0.97 

0.85 

102.3 

6.3 

0.12 

89.9 

100.6 

0.0 

4.5 

0.00 

1.5 

Grid 

0.00 

0.01 

100.1 

7.2 

0.08 

100.0 

100.1 

20.0 

2.2 

0.07 

0.05 

Lattice 

0.00 

0.02 

100.0 

5.6 

0.17 

100.0 

100.0 

20.0 

1.9 

0.13 

0.02 

Noise-added 

Lines* 

0.01 

0.48 

103.0 

7.4 

0.08 

97.5 

96.5 

12.6 

4.3 

0.01 

2.8 

Lines** 

0.71 

0.69 

101.2 
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0.20 

100.0 

98.6 

13.7 
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0.10 

3.6 

Grid 

-0.36 

-0.01 

101.0 

7.7 

0.07 

100.1 

98.2 

15.5 

2.1 

0.07 

2.9 

Lattice 

0.81 

-0.47 

99.6 

5.6 

0.16 

95.3 

101.1 

18.6 

1.3 

0.39 

2.2 

*  lines  pass  directly  over  both  of  the  targets,  required  tighter  limits  on  XY  to  resolve  small  target 
**  lines  pass  10cm  offset  from  center  of  both  targets 
susceptibility  conversion:  k%\  =  4n  kcgs 


Figure  6:  Screen  capture  from  the  Encom  ModelVision  software.  It  includes  2D  plan,  2D 
profile,  and  3D  data-surface  visualization  capabilities  and  inversion  parameter  control. 


Conclusion 

Three-dimensional  processing  and  presentation  of  geophysical  data  collected  with  3D 
positioning  offers  several  advantages.  The  primary  benefit  is  the  ability  to  improve  noise  levels 
by  spatially  correlating  data  in  three  dimensions  without  forcing  the  data  onto  a  common  plane. 
A  new  visualization  technique  using  gradient-strings  was  introduced  to  enable  a  full  3D  lattice  to 
be  reduced  to  its  most  intrinsic  qualities  for  interpretation.  For  a  string  to  be  a  representation  of  a 
true  ground  source,  it  must  (a)  exceed  an  intensity  threshold  (b)  be  continuous  over  several  layers 
and  (c)  increase  in  intensity  with  depth.  By  eliminating  strings  that  do  not  meet  these  criteria, 
the  result  is  a  simple  line  through  the  3D  data  space  that  points  directly  to  the  target(s).  To 
further  aid  in  interpretation,  the  line  is  color-coded  to  gradient  intensity. 

The  effects  of  second  order  gradients  were  also  briefly  demonstrated.  The  response  of 
the  larger  target  was  visible  at  all  heights,  but  the  gradient  of  the  shallower  target  increased  much 
more  quickly  with  proximity.  Since  cued  investigations  typically  use  analog  instruments  held 
very  close  to  the  ground,  it  is  a  debatable  and  subjective  call  as  to  which  target  would  actually  be 
excavated  by  a  field  crew  -  the  larger  target,  or  the  stronger  response.  This  digital  3D  approach 
would  make  it  much  simpler  for  field  crews  to  assess  multiple  targets  in  close  proximity  to  the 
single  dig  location  they  may  be  provided. 


Inversion  results  of  3D  data  were  shown  to  be  equivalent  to  an  equal  density  of  highly 
controlled  2D  data,  but  significantly  better  than  typical  line  survey  results.  The  combination  of 
gradient  strings  for  simple  visualization  and  the  numerical  inversion  provide  mutual  feedback 
into  the  decision  making  process.  The  visualization  and  manipulation  of  gradient  strings  will 
allow  the  operator  to  refine  processing  parameters  and  to  deselect  targets  from  noise  or  to  pin¬ 
point  multiple  targets  within  close  proximity  to  one  another.  In  return,  the  inversion  will  provide 
quantitative  results  in  the  form  of  depths,  sizes  and  fit  coefficients. 

Acknowledgements 

Funding  for  this  work  was  provided  as  a  SEED  project  by  the  SERDP  UXO  office  under 
Anne  Andrews.  Oak  Ridge  National  Laboratory  is  managed  by  UT-Battelle,  LLC  for  the  U.S. 
Department  of  Energy  under  contract  DE-AC05-00OR22725.  The  submitted  manuscript  has 
been  authored  by  a  contractor  of  the  U.S.  Government.  Accordingly,  the  U.S.  Government 
retains  a  nonexclusive,  royalty-free  license  to  publish  or  reproduce  the  published  fonn  of  this 
contribution,  or  allow  others  to  do  so  for  U.  S.  Government  purposes. 


Appendix  A4 

Oral  Presentation  -  SAGEEP  2004 
February  2004 


'-^3D  Processing  and 
Visualization  of  3D  Magnetic 

Data 


T.  Jeffrey  Gamey,  William  E.  Doll,  Abraham  Emond,  Les  P.  Beard 
Oak  Ridge  National  Laboratory,  Oak  Ridge,  TN 


Background 

>  OE  Techs  don’t  trust  digital  geophysics 

>  They  DO  trust  analog  geophysks  ! ! 

>  Decision  making  is  taken  out  of  their  hands 

>  Ultimate  objective  is  to  simplify  the 
visualization  and  discrimination  tools  to  the 
point  where  Techs  can  use  and  trust  them  to 
make  decisions  that  will  speed  clearance 

>  Still  retains  digital  record  of  each  decision 
as  required  by  regulators 


1 


'Overview 


>  Definitions 

\ 

>  Traditional  Methods 

>  State-of-the-Art  Methods 

>  3D  Method 

>  Generalized  example  with  modeled  data 

>  Gradient  strings 

>  3D  vs  2D  inversion  comparison 

>  Conclusions 


Collecti 


■^Definitions: 
^Presentation  -  Inversion 


>  Collection  refers  to  number  of  dimension  coordinates 
collected  with  the  geophysical  data 

•  rarely  ID  (distance  down  line) 

•  usually  2D  (includes  XY) 

•  recently  3D  (XYZ) 

>  Z  measurement  techniques  include  DGPS,  laser,  radar, 
ultra-sonic 

>  Z  accuracy  varies  with  technology;  DGPS  Z  measure  is 
roughly  half  as  good  as  XY 

>  Z  reference  planes  vary  with  technology;  HAE  is  good  for 
inversion,  not  so  good  for  ground  follow-up 
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~ — ——JDefinitions: 

Collection  -  -  Inversion 

>  Presentation:  number  of  dimension  axes  on  map, 
profiles  (ID),  grids  (2D),  solids  (3D) 

•  same  definition  applies  to  processing  techniques 

•  note  some  “3D  maps”  use  geophysical  axis  for  Z 

•  note  limitations  of  2D  media  (paper,  screens) 

>  Presenting  third  axis  (dimension  or  geophysical)  on 
2D  media  uses  perspective  and  shading  to  create 
illusion  of  three  dimensions 

>  Presenting  “4D”  data  requires  additional  techniques 
to  visualize  the  interior  of  the  data  space 

(slice,  fence,  transparency,  iso-surface) 


' — ——-Definitions: 

Collection!  Pre&er|ation  - 

>  Inversion  refers  to  number  of  variable  dimension 
coordinates  in  geophysical  data  used  in  target 
calculation 

>  Standard  techniques  to  locate  discrete  targets  such 
as  Analytic  Signal  peaks  and  Euler  are  based  on 
grids  with  2D  positions  (fixed  Z,  very  fast,  simple) 

>  Newer  techniques  use  a  collection  of  points  (line  or 
set  of  lines)  with  measured  Z  positions  (3D) 
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Traditional  Methods 

>  Data  are  collected  in  2D  at  nominal  height  AGL 

>  Then  processed  in  ID  (freq  filters), Ridded, 
additional  processing  in  2D  (filters,  AS) 

>  Then  interpreted  in  2D  (peaks,  Euler) 

:  Fast-simple-intuitive,  relates  anomalies  on 
adjacent  lines  for  discrimination  of  target  shapes 

:  Ignores  differences  in  collection  height;  forces 
all  data  onto  common  plane;  converts  height 
changes  into  grid  artifacts;  spends  more  effort 
removing  height  artifacts  than  improving  quality 
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State-of-the-Art  Methods 

>  Data  are  collected  in  3D  at  nominal  height  AGL 

>  Then  processed  in  ID  (freq  filters),  gridded, 
additional  processing  in  2D  (filters,  AS;  peaks) 

>  Then  interpreted  in  3D  from  pre-grid  data 

:  Uses  height  data,  bypass  grid  artifacts 

:  Misses  opportunity  to  reduce  noise  for  data 
inversion  through  cross-line  correlation;  target 
selection  still  distorted  by  forcing  height  differences 
onto  common  plane 


State-df-the-Art  Methods 

>  Data  collected  in  3D  at  nominal  height  AGL 

>  Gridded  data  still  show  artifacts,  bqt  inversion 
avoids  these  problems  by  stepping  back  to  use 
profile  data  only 
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3D  Method 

>  Data  are  collected  in  3D  at  random  heights  AGL 

>  Then  processed  in  ID  (filters?),  gridded  (3D), 
additional  processing  in  3D  (filters,  A5>)  strings) 

>  Then  interpreted  in  3D  from  3D  gridded  data 

:  Uses  height  data  and  three  dimensional 
correlation  to  reduce  noise  prior  to  target  selection 
and  inversion 

:  Practical  only  for  cued  investigations; 
integrated  system  not  yet  available 


Example  of  3D  Method 

>  Modeled  response  from  two  small  spherical 
targets  with  different  radii  and  burial  depths, 
simulating  60mm  mortar  and  20mm  shell 

>  Data  were  collected  at  random  positions  within 
3D  space  (100cm  x  100cm  x  50cm) 

>  Sampled  at  l/6th  full  density  (same  number  of  data 
points  as  fully  populated  grid  at  10cm  interval) 

>  Sample  heights  from  10-60cm  AGL 

>  Random  noise  added  to  positions  and  magnetic 
data  for  testing  purposes 
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Example  of  3D  Method 


>  Data  are  collected  in  3D  at  random  heights  AGL 


Top 


Base 


Example  of  3D  Method 

>  Data  are  “gridded”  in  3D  onto  a  lattice  of  regularly 
spaced  points,  creating  a  solid  volume  of  data 

>  Six  different  3D  gridding  techniques  were  tested 
using  Rockware  SolidModel  routines 

>  “anisotropic  inverse-distance  weighting”  produced 
the  best  results  based  on  least  error  between 
interpolated  and  modeled  response  regardless  of 
degree  of  interpolation  tested 

>  3D  Hanning  (cosine  roll-off)  and  triangular  (linear 
roll-off)  filters  were  designed  to  supplement  square 
(vertical  roll-off)  filter  provided 
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Example  of  3D  Method 

>  Data  are  “gridded”  in  3D  ontoa  lattice  of  regularly 
spaced  points,  creating  a  solid  volume  of  data 


Example  of  3D  Method 


>  Visualization  not  required  for  inversion,  but  it  is 
required  for  interpretation 

>  Interpreter  must  at  least  specify  initial  target  shape 
(sphere,  ellipsoid,  continuous  pipeline)  and  whether  to 
invert  for  single  or  multiple  targets 

>  Additional  visualization  techniques  are  required  for  3D 

•  Solids  only  show  surface  readings 

•  Total  field  after  3D  Hanning 
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Example  of  3D  Method 


>  Additional  visualization  techniques  are  required  for  3D 

•  grids 

•  slices 

•  fences 


Example  of  3D  Method 


>  Additional  visualization  techniques  are  required  for  3D 

•  iso-surfaces 
(<50,000nT 
and 

>50,002nT) 
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Gradient  Strings 

>  Isolate  individual  peaks  and  display 

>  Numerous  peaks  will  be  caused  by  noise 


Top  Top 


all  strings  strings  >  2pts 
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Gradient  Strings 

>  Number  of  strings  can  be  reduced  by  length, 
strength  and  continuity 


Top  Top 


strings  >  0.3nT/m 


strings  >  2pts 
and  >  0.3nT/m 


Inversion  Comparison 

>  Modeled  data  from  two  targets  based  on  traditional 
survey  lines  and  3D  data  collection 

>  Survey  lines  at  0.5m  spacing  0.4m  height;  one  set 
directly  over  targets,  one  offset  by  10cm 

>  Line  data  required  5:1  interpolation  to  fill  2D  grid, 
but  still  twice  the  density  of  typical  lm  survey 

>  3D  data  collected  at  random  locations 

>  3D  data  required  6:1  interpolation  to  fill  3D  solid 
(same  number  of  data  points  as  2D  grid) 

>  Noise-free  and  noise-added  sets  compared 


I n ve rsion  Co mparison 

>  Line  data  detection  of  smaller  target  problematic 
(insufficient  density) 

>  Inversion  of  line  data  had  trouble  resolying  location  of 
smaller  target,  especially  if  line  offset  from  target 

>  3D  strings  detected  both  targets 

>  Inversion  of  3D  data  showed  better  fit  results 

(line  data  2cm  off  large  target,  3D  exact  fit  on  both) 

>  Adding  positioning  noise  evened  up  the  comparison  by 
degrading  quality  of  3D  faster  than  line  data  due  to 
higher  gradients  at  lower  heights,  but  3D  still  better 
(line  data  3/10cm  off,  3D  l/5cm  off) 

>  3D  resolved  target  size  better  (10%  vs  50-150%) 


Conclusions 


>  Traditional  collection  methods  suffer  from  gridding 
artifacts  caused  by  height  variations 

>  Modern  inversion  methods  eliminate  this  by  collecting 
3D  positions  and  inverting  collections  of  lines,  but  this 
loses  the  noise  reduction  capabilities  of  cross-line 
correlation 

>  3D  data  collection  solves  both  of  these  problems  plus 
it  improves  inversion  results,  but  is  practical  only  for 
cued  investigations 

>  Gradient  strings  were  introduced  to  aid  visualization 
of  3D  solids  prior  to  setting  inversion  parameters 

>  The  technique  is  as  applicable  to  EM  as  to  magnetic  I 
methods 
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3D  Method 

>  Data  are  collected  in  3D  at  random  heights  AGL 
then  gridded  in  3D  for  processing  and  interpretation 
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Appendix  B1 

Color  Contours  with  Targets 


Color  contours  of  2D  Total  Field  Magnetic  data  with  synthetic  targets  in  their  true 
locations  and  orientation.  Horizontal  units  in  cm,  total  field  units  in  nT. 
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Vlodel  3  (spheres  at  depths  0.6m  and  0.05m) 
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Model  5  (Horizontal  ellipsoid  at  0.6m) 
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Model  7  (Horizontal  ellipsoids  at  depths  0.6m  and  0.05m) 


VI odd  8  (Vertical  ellipsoids  at  depths  0.6m  and  0.05m) 


Vlodel  9  (Randomly  oriented  ellipsoids  at  depth  0.6m  and  0.2m) 


Appendix  B2 
Gradient  Strings 

Gradient  Strings  of  3D  Analytic  Signal  data.  Data  space  is  2m  x  2m  x  0.5m.  Viewed 
above  from  the  SE.  Analytic  Signal  units  in  nT/m. 
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Model  4  (Horizontal  ellipsoid  at  depth  0.2m) 
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Model  6  (Vertical  ellipsoid  at  depth  0.6m) 
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Model  7  (Horizontal  ellipsoids  at  depths  0.6m  and  0.05m) 


< 


Model  8  (Vertical  ellipsoids  at  depths  0.6m  and  0.05m) 
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Model  9  (Randomly  oriented  ellipsoids  at  depth  0.6m  and  0.2m) 


Appendix  B3 

UXAnalyze  Dipole  Inversion  Results 

This  code  is  currently  under  development  by  AETC  and  Geosoft  under  a  concurrent 
ESTCP  project.  The  authors  would  like  to  extend  their  gratitude  to  AETC  for  running 
these  models  for  us  when  our  original  attempts  to  run  the  older  IDL-based  MTADS-DAS 
code  failed  to  invert  several  of  the  2D  input  data  sets.  The  positional  accuracy  of  the  two 
runs  are  comparable,  but  the  UXAnalyze  version  inverted  the  2D  data  with  greater 
success. 

As  with  the  IDL-based  code,  this  code  operates  on  one  anomaly  at  a  time.  Where 
signatures  overlap,  data  must  be  manually  divided  between  anomalies  based  on  position. 
For  the  data  in  this  Appendix,  no  segmentation  of  anomalies  was  attempted.  Targets 
generally  appear  as  single  asymmetric  anomalies  in  the  2D  data,  and  the  3D  data  cannot 
be  satisfactorily  gridded.  As  a  result,  all  Models  are  treated  as  single  targets  for  purposes 
of  inversion  and  results  should  be  compared  to  the  “combined”  targets  presented  in  the 
body  of  the  text.  Position  errors  are  reported  as  the  linear  distance  from  the  detection 
point  to  the  closest  true  target  location  including  both  horizontal  and  vertical  components. 

This  inversion  process  is  essentially  the  same  as  the  “combined”  target  inversion 
presented  in  the  body  of  the  text.  When  these  “combined”  targets  are  taken  in  to 
consideration,  the  results  both  in  terms  of  fit  and  accuracy  are  virtually  the  same  for  both 
the  IDL  and  UXAnalyze  codes. 


Table  B3  - 1 :  2D  data 

It  coefficient. 

Fit  (2D  data) 

Fit  (3D  data) 

Model  1 

0.972 

0.000 

Model  2 

0.928 

0.925 

Model  3 

0.939 

0.497 

Model  3t 

0.960 

0.002 

Model  4 

0.957 

0.835 

Model  5 

Not  strong  enough 

0.495 

Model  6 

0.983 

Too  much  change 

Model  7 

Not  strong  enough 

0.331 

Model  8 

0.952 

0.628 

Model  9 

0.975 

0.570 

Table  B3-2:  2D  and  3D  dipole  inversion  results.  Average  location  error  of  0.16m  with 
two  failures  for  2D  input  data.  Average  location  error  of  0.29m  with  one  failure  for  3D 
input  data. _ _ _ _ 


2D  data 

Error 

3D  data 

Error 

X 

Y 

Z 

size 

X 

Y 

Z 

size 

True 

-0.10 

-0.20 

0.60* 

0.10 

-0.10 

-0.20 

0.60* 

0.10 

Model  1 

-0.17 

-0.21 

0.15 

0.05 

0.09 

-0.05 

-0.28 

-0.12 

0.00 

0.33 

Model  2 

-0.03 

-0.17 

0.52 

0.04 

0.11 

-0.16 

-0.26 

0.43 

0.04 

0.19 

Model  3 

-0.06 

-0.02 

0.50 

0.05 

0.21 

0.00 

-0.26 

0.78 

0.05 

0.21 

Model  3t 

-0.08 

-0.02 

0.53 

0.05 

0.19 

0.10 

0.36 

-0.15 

0.00 

0.21  J 

Model  4 

-0.07 

-0.18 

0.21 

0.05 

0.04 

-0.17 

-0.24 

0.14 

0.04 

0.10 

Model  5 

— 

— 

-0.33 

-0.24 

0.19 

0.02 

0.47 

Model  6 

-0.06 

-0.11 

0.49 

0.06 

0.15 

— 

— 

Model  7 

— 

— 

0.02 

0.35 

0.00 

0.01 

o.i  ij 

Model  8 

-0.02 

0.03 

0.49 

0.06 

0.27 

0.00 

0.17 

0.83 

0.07 

0.45 

Model  9 

-0.01 

0.03 

0.48 

0.06 

0.28 

0.02 

0.23 

0.51 

0.05 

0.45 

*depth  of  models  1  and  4  is  0.20m 

^position  error  with  respect  to  20mm  target  at  (0.10,  0.30,  0.05). 
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Measured  Data 

Fit  Results 
X  (m):  7.99 

Y  (m):  12.03 

Depth  (m):  0.48 

Inclination  (°):  56.1 
Declinaton  (°):  50.5 
Size  (m):  0.057 

Coherence:  0.975 

Solid  Anglef):  25.7 


Modeled  Data 


2D  -  Anomaly  10 


Figure  B3-1 :  Screen  capture  of  UXAnalyze  results  from  Model  9. 


Appendix  B4 

Solid-body  Inversion  of  Interpolated  Data 


In  addition  to  the  solid-body  inversion  of  the  2D  lines  and  3D  “cloud”  of  data  points, 
inversion  was  also  run  on  the  interpolated  2D  grids  and  3D  lattices.  The  results  are 
significantly  skewed  by  the  introduction  of  interpolation  errors.  This  impacted  the  3D 
data  more  than  the  2D  data  due  to  the  poorer  algorithms  and  the  greater  requirement  for 
interpolation.  The  2D  data  required  interpolation  from  205  to  441  points  (1:2),  whereas 
the  3D  data  required  interpolation  from  300  to2,205  points  (1:7). 

The  weakness  in  the  Rockware  gridding  algorithm  (used  for  all  3D  data)  is  demonstrated 
using  Model  9rw.  In  this  case,  the  Rockware  routines  were  used  to  grid  the  2D  data,  and 
can  be  compared  directly  to  the  results  using  the  Geosoft  routines  as  shown  in  the  2D 
Model  9  results.  The  entire  2D  data  set  was  not  treated  with  the  Rockware  processing  in 
order  to  illustrate  the  perfonnance  of  the  3D  method  against  the  current  standard  of 
processing  algorithms. 

The  Model  9t  reported  in  the  3D  results  is  the  fully-populated  noise-free  version  of 
Model  9.  Since  this  model  is  fully-populated  it  does  not  rely  on  interpolation  of  any  kind. 
The  errors  in  this  case  are  strictly  the  result  of  the  mis-match  in  treatment  of  ellipsoids 
between  the  MAGMOD  forward  modeling  code  and  the  ModelVision  inversion  code. 

Note  that,  unlike  the  tables  in  the  body  of  the  report,  the  positioning  errors  reported  here 
are  the  z-error  and  not  the  radial-error. 


Table  B4-1 :  2D  RMS  error  w 

len  fitting  a  sphere  or  ellipsoic 

to  the  data. 

RMS  error  (sph) 

RMS  error  (ellip) 

Model  1 

2.30 

2.21 

Model  2 

1.80 

1.63 

Model  3 

1.82 

1.64 

Model  3t 

0.01 

0.08 

Model  4 

5.70 

5.68 

Model  5 

8.76 

8.67 

Model  6 

2.58 

2.56 

Model  7 

8.32 

8.94 

Model  8 

2.40 

2.36 

Model  9 

4.71 

4.60 

Model  9rw 

6.12 

5.26 

Table  B4-2:  2D  Spherical  inversion  results  (fitting  a  spherical  model) 


60mm 

z-error 

20mm 

z-error 

X 

Y 

Z 

rad 

X  Y 

Z 

rad 

True 

-0.10 

-0.20 

0.60* 

0.10 

0.10  0.30 

0.05 

0.03 

Model  1 

-0.05 

-0.19 

0.24 

0.05 

0.04 

Not  modeled 

— 

Model  2 

-0.07 

-0.20 

0.63 

0.05 

0.03 

Not  modeled 

— 

Model  3 

-0.07 

-0.21 

0.64 

0.05 

0.04 

-0.02  0.29 

0.15 

0.02 

0.10 

Model  3t 

-0.10 

-0.20 

0.60 

0.05 

0.00 

0.10  0.30 

0.05 

0.02 

0.00 

Model  4 

-0.04 

-0.38 

0.23 

0.04 

0.03 

Not  modeled 

— 

Model  5 

-0.05 

-0.50 

0.56 

0.04 

0.04 

Not  modeled 

— 

Model  6 

-0.08 

-0.10 

0.65 

0.06 

0.05 

Not  modeled 

— 

Model  7 

-0.05 

-0.50 

0.60 

0.04 

0.00 

-0.01  0.10 

0.03 

0.01 

0.02 

Model  8 

-0.07 

-0.14 

0.62 

0.06 

0.02 

-0.03  0.35 

0.32 

0.03 

0.27 

Model  9 

-0.20 

-0.12 

0.62 

0.05 

0.02 

-0.09  0.29 

0.34 

0.03 

0.29 

Model  9rw 

-0.23 

-0.20 

0.50 

0.05 

0.10 

-0.10  0.23 

0.43 

0.04 

0.38 

*depth  of  models  1  and  4  is  0.20m 


Table  B4-3:  2D  Ellipsoid  inversion  results  (positions) 


60mm 

z-error 

20mm 

z-error 

X 

Y 

Z 

X  Y 

Z 

True 

-0.10 

-0.20 

0.60* 

0.10  0.30 

0.05 

Model  1 

-0.06 

-0.20 

0.25 

0.05 

Not  modeled 

— 

Model  2 

-0.10 

-0.30 

0.65 

0.05 

Not  modeled 

— 

Model  3 

-0.08 

-0.21 

0.63 

0.03 

-0.02  0.28 

0.16 

0.11 

Model  3t 

-0.10 

-0.20 

0.61 

0.01 

0.10  0.30 

0.05 

0.00 

Model  4 

-0.04 

-0.38 

0.23 

0.03 

Not  modeled 

— 

Model  5 

-0.04 

-0.50 

0.56 

0.04 

Not  modeled 

— 

Model  6 

-0.08 

-0.10 

0.65 

0.05 

Not  modeled 

— 

Model  7 

-0.04 

-0.43 

0.57 

0.03 

-0.10  0.10 

2.00 

1.95 

Model  8 

-0.08 

-0.16 

0.61 

0.01 

-0.03  0.34 

0.36 

0.21 

Model  9 

-0.22 

-0.10 

0.63 

0.03 

-0.06  0.27 

0.31 

0.14 

Model  9rw 

-0.27 

-0.12 

0.68 

0.08 

-0.04  0.19 

0.26 

0.34 

*depth  of  models  1  and  4  is  0.20m 


Table  B4-4:  2D  Ellipsoid  inversion  results  (dimensions) 


60mm 

error 

20mm 

Error 

A 

B 

C 

%Vol 

A  B 

C 

%Vol 

True 

0.051 

0.051 

0.051 

0.017  0.017 

0.017 

Model  1 

0.12 

0.03 

0.03 

20 

Not  modeled 

— 

Model  2 

0.28 

0.02 

0.02 

25 

Not  modeled 

— 

Model  3 

0.26 

0.02 

0.02 

29 

0.04  0.01 

0.01 

30 

Model  3t 

0.10 

0.03 

0.03 

29 

0.02  0.02 

0.02 

28 

True 

0.100 

0.030 

0.030 

0.035  0.010 

0.010 

Model  4 

0.09 

0.03 

0.03 

21 

Not  modeled 

— 

Model  5 

0.17 

0.03 

0.01 

33 

Not  modeled 

— 

Model  6 

0.57 

0.15 

0.01 

162 

Not  modeled 

— 

Model  7 

1.00 

0.03 

0.01 

27 

0.02  0.02 

0.01 

71 

Model  8 

0.19 

0.06 

0.05 

104 

0.12  0.03 

0.02 

1145 

Model  9 

0.18 

0.06 

0.05 

83 

0.04  0.03 

0.01 

701 

Model  9rw 

0.13 

0.06 

0.03 

158 

0.04  0.04 

0.02 

841 

Table  B4-5:  3D  RMS  error  w 

len  fitting  a  sphere(s)  or  ellipsoid(s)  to  the  data. 

Error  (sph) 

Error  (ellip) 

Model  1 

1.79 

1.68 

Model  2 

3.57 

3.49 

Model  3 

1.91 

2.36 

Model  3t 

0.01 

0.01 

Model  4 

2.43 

2.36 

Model  5 

4.99 

4.92 

Model  6 

3.17 

3.55 

Model  7 

3.44 

1.99 

Model  8 

1.98 

2.35 

Model  9 

1.73 

1.72 

Model  9t 

0.66 

0.63 

Table  B4-6:  3D  Spherical  inversion  results  (fitting  a  spherical  model) 


60mm 

z-error 

20mm 

z-error 

X 

Y 

z 

rad 

X 

Y 

Z 

rad 

True 

-0.10 

-0.20 

0.60* 

0.10 

0.10 

0.30 

0.05 

0.03 

Model  1 

-0.11 

-0.24 

0.30 

0.06 

0.10 

Not  modeled 

— 

Model  2 

-0.14 

-0.30 

0.86 

0.06 

0.26 

Not  modeled 

— 

Model  3 

-0.14 

-0.30 

0.90 

0.07 

0.30 

0.08 

0.27 

0.13 

0.02 

0.08 

Model  3t 

-0.10 

-0.20 

0.60 

0.05 

0.00 

0.10 

0.30 

0.05 

0.02 

0.00 

Model  4 

-0.12 

-0.39 

0.35 

0.06 

0.15 

Not  modeled 

— 

Model  5 

-0.15 

-0.50 

1.10 

0.07 

0.50 

Not  modeled 

— 

Model  6 

-0.13 

-0.20 

0.76 

0.08 

0.16 

Not  modeled 

— 

Model  7 

-0.15 

-0.50 

0.78 

0.05 

0.18 

0.22 

0.21 

0.14 

0.01 

0.09 

Model  8 

-0.14 

-0.22 

0.77 

0.08 

0.17 

0.09 

0.30 

0.15 

0.03 

0.10 

Model  9 

-0.30 

-0.16 

0.87 

0.08 

0.27 

0.04 

0.24 

0.11 

0.02 

0.06 

Model  9t 

-0.20 

-0.10 

0.58 

0.06 

0.02 

0.07 

0.28 

0.05 

0.02 

0.00 

*depth  of  models  1  and  4  is  0.20m 


Table  B4-7:  3D  Ellipsoid  inversion  results  (positions) 


60mm 

z-error 

20mm 

z-error 

X 

Y 

Z 

X  Y 

Z 

True 

-0.10 

-0.20 

0.60* 

0.10  0.30 

0.05 

Model  1 

-0.11 

-0.25 

0.25 

0.05 

Not  modeled 

— 

Model  2 

-0.15 

-0.31 

0.87 

0.27 

Not  modeled 

— 

Model  3 

-0.12 

-0.24 

0.92 

0.32 

0.08  0.29 

0.02 

0.03 

Model  3t 

-0.10 

-0.20 

0.60 

0.00 

0.10  0.30 

0.05 

0.00 

Model  4 

-0.12 

-0.41 

0.36 

0.16 

Not  modeled 

— 

Model  5 

-0.18 

-0.50 

1.09 

0.49 

Not  modeled 

— 

Model  6 

-0.14 

-0.22 

1.10 

0.50 

Not  modeled 

— 

Model  7 

-0.27 

-0.50 

0.91 

0.31 

0.07  0.22 

0.07 

0.02 

Model  8 

-0.17 

-0.31 

0.86 

0.26 

0.03  0.22 

0.35 

0.30 

Model  9 

-0.30 

-0.16 

0.88 

0.28 

0.04  0.25 

0.10 

0.05 

Model  9t 

-0.22 

-0.10 

0.60 

0.00 

0.06  0.29 

0.05 

0.00 

*depth  of  models  1  and  4  is  0.20m 


Table  B4-8:  3D  Ellipsoid  inversion  results  (dimensions) 


60mm 

error 

20mm 

Error 

A 

B 

C 

%Vol 

A  B 

C 

%Vol 

True 

0.051 

0.051 

0.051 

0.017  0.017 

0.017 

Model  1 

0.29 

0.07 

0.01 

50 

Not  modeled 

— 

Model  2 

0.24 

0.04 

0.03 

102 

Not  modeled 

— 

Model  3 

0.12 

0.06 

0.05 

147 

0.01  0.01 

0.01 

80 

Model  3t 

0.07 

0.07 

0.03 

1 

0.02  0.02 

0.02 

0 

True 

0.100 

0.030 

0.030 

0.035  0.010 

0.010 

Model  4 

0.44 

0.04 

0.04 

122 

Not  modeled 

— 

Model  5 

0.25 

0.14 

0.01 

306 

Not  modeled 

— 

Model  6 

0.57 

0.15 

0.01 

872 

Not  modeled 

— 

Model  7 

1.00 

0.03 

0.01 

258 

0.02  0.02 

0.01 

21 

Model  8 

0.19 

0.06 

0.05 

433 

0.12  0.03 

0.02 

1659 

Model  9 

0.18 

0.06 

0.05 

503 

0.04  0.03 

0.01 

208 

Model  9t 

0.20 

0.08 

0.02 

179 

0.06  0.01 

0.01 

153 

The  spherical  inversion  of  the  spherical  targets  slightly  out-performed  the  other 
combinations  (e.g.  spherical  inversion  of  ellipsoidal  targets).  The  most  disturbing  result 
is  the  very  poor  (x,  y)  positioning  the  noise-free  ellipsoid  data  (Model  9t),  even  though 
the  depth  is  accurately  calculated. 

Inversion  using  the  2D  Geosoft  grids  showed  a  slightly  poorer  overall  fit,  especially  for 
Models  5  and  7.  The  ellipsoid  inversion  for  the  20mm  in  Model  7  should  be  classified  as 
a  non-detection  since  the  results  pushed  the  target  to  the  extreme  limits  set  for  distance 
and  depth  (Table  B4-3,  blue  text).  The  radial  distance  to  the  target  center  was 
comparable  between  the  2D  and  3D  data  sets,  averaging  0.31  and  0.28m  respectively,  but 
the  depth  error  in  the  2D  data  was  nearly  an  order  of  magnitude  better  than  the  3D  data. 
The  size  estimates  based  on  target  volume  were  also  about  twice  as  good  for  the  2D  data. 
This  is  attributed  to  the  superiority  of  the  2D  Geosoft  “minimum  curvature”  gridding 
algorithms  over  the  3D  “inverse-distance”  algorithm  in  Rockware.  This  is  supported  by 
the  considerably  poorer  perfonnance  of  the  Model  9rw.  In  this  case,  the  2D  data  were 
gridded  and  processed  using  the  Rockware  code  rather  than  the  standard  Geosoft  code. 


Appendix  B5 

Rockware  Gridding  Algorithms 


These  descriptions  are  taken  directly  from  the  Rockware  help  manual  provided  by  the 
software  company.  They  cover  the  six  gridding  algorithms  tested  under  this  project. 


1)  Closest-Point  gridding  algorithm 

Click  in  this  radio  button  to  select  the  Closest-Point  algorithm,  in  which  the  value  of  a 
voxel  node  is  set  to  be  equal  to  the  value  of  the  nearest  data  point,  regardless  of  its 
distance  from  the  point  or  the  value  of  its  other  neighbors. 

This  method  is  useful  when  generating  models  in  which  the  values  are  not  gradational. 


2)  Isotropic  Distance  gridding  algorithm 

Click  in  the  Inverse-Distance/Isotropic  radio  button  to  select  the  Inverse-Distance 
algorithm  which  assigns  a  voxel  node  value  based  on  the  weighted  average  of 
neighboring  data  points.  Unlike  the  "Anisotropic"  method,  this  method  uses  all  data 
points  in  interpolating  node  values. 

During  the  interpolation  of  solid  model  nodes,  the  value  of  each  data  point  is  weighted 
according  the  inverse  of  its  distance  from  the  node,  taken  to  the  second  power  (Inverse- 
Distance  squared).  By  raising  this  to  an  exponent  of  “2”,  this  method  assures  that  distant 
points  will  have  less  influence  on  the  node  value  than  closer  points. 

This  method  is  useful  when  modeling  uniformly  distributed  data  in  non-stratiform 
environments. 


3)  Anisotropic  Distance  gridding  algorithm 

Click  in  this  radio  button  to  select  the  Octant-Based  Inverse-Distance  algorithm,  which 
assigns  a  voxel  node  value  based  on  the  weighted  average  of  neighboring  data  points. 
Unlike  the  "Isotropic"  method,  this  method  includes  only  the  closest  data  point  within 
each  octant  around  the  node  when  computing  the  inverse-distance  weighted  average. 

This  is  useful  for  modeling  drill-hole  based  data  in  stratiform  deposits. 


4)  Weighted  Distance  gridding  algorithm 

Click  in  this  radio  button  to  select  the  Inverse-Distance  algorithm  that  will  use  all  of  the 
available  data  points  when  computing  a  node’s  value. 

But,  instead  of  automatically  using  a  weighting  exponent  of  “2”,  the  program  allows  the 
user  to  assign  different  weighting  exponents  to  control  points  oriented  vertically  versus 
horizontally  from  the  node.  The  greater  the  exponent  you  enter,  the  less  influence  those 
data  points  will  have. 

Expand  this  heading  to  establish  the  weighting  values. 


5)  Directional  Weighted  gridding  algorithm 

Click  in  the  Directional  Weighting  radio  button  to  select  the  Directional- Weighting 
algorithm.  This  method  will  bias  the  modeling  in  a  user-specified  direction. 

Expand  this  heading  to  specify  the  direction  and  strength. 


6)  Vertical  Biasing  gridding  algorithm 

This  modeling  method  uses  Inverse-Distance  modeling,  in  which  the  value  of  each  of  the 
data  points  is  weighted  according  to  the  inverse  of  its  distance  from  the  voxel  node  being 
computed,  taken  to  a  user-selected  power.  The  difference  with  Horizontal  Biasing  is  that 
the  user  can  define  a  vertical  distance  from  each  voxel  node  beyond  which  points  will  no 
longer  be  used  in  computing  the  node  value.  This  allows  the  user  to  horizontally  bias  the 
interpolation  process  such  that  stratiform  deposits  such  as  sand  and  gravel  may  be 
modeled  more  effectively 


Expand  this  heading  to  establish  the  biasing  "bubble. 


